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ABSTRACT 

The  solution  of  a  beam  on  an  elastic  nonlinear  foundation  by  the 
Finite  Element  Method  is  presented  in  this  thesis.     Discontinuous 
(Winkler)  and  continuous  foundations  are  considered.      The  general 
formulation  is  based  on  the  Galerkin  method.      The  solution  technique 
for  the  linear  case  uses  Gauss  elimination  and  the  solution  technique 
for  the  nonlinear  case  is  based  on  Brown's  method  (a  modified  Newton- 
Raphson).     Some  illustrative  examples  are  presented.     A  brief  com- 
parison of  continuous  and  non- continuous  foundations  is  made. 
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I.     INTRODUCTION 

In  treating  the  various  problems  related  to  structures  on  elastic 
foundations,    it  has  been  customary  for  engineers  to  assume  that  the 
pressure  in  the  foundation  is  linearly  proportional  at  every  point  to  the 
deflection  of  the  foundation  at  that  point,    independent  of  the  pressure  or 
deflection  occurring  in  other  parts  of  the  foundation.      This  assumption, 
is  mathematically  by  far  the  simplest  that  one  can  make  regarding  the 
nature  of  a  supporting  elastic  medium.     It  assumes  a  complete  lack  of 
continuity  in  the  material  of  the  foundation  as  if  it  consisted  of  a  series 
of  independent  springs  which  deflect  when  directly  loaded.      This  theory 
has  proved  adequate  for  calculating  stresses  and  deflections  in  railroad 
tracks  and  has  found  many  applications  to  plate  and  shell  structures 
[Ref.    1].     However  no  particular  claim  has  been  made  that  the  deforma- 
tion or  pressure  distribution  in  actual  earth  foundations  could  be  pre- 
dicted by  this  method. 

In  recent  years,    some  authors  have  treated  foundations  as  a  con- 
tinuous medium  which,    in  contrast  to  the  Winkler  hypothesis,    repre- 
sents the  case  of  complete  continuity  in  the  material  [Ref.    9].     It  is 
evident  that  the  two  types  of  foundation  naodels  lead  to  different  results, 
but  how  quantitatively  significant  is  the  difference  has  not  previously 
been  made  clear. 
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The  properties  of  actual  foundations,    however,    seem  to  lie  some- 
where in  the  gap  between  these  two  extreme  cases.      The  physical  prop- 
erties of  soils  are  obviously  of  a  very  complicated  nature.     Their  de- 
formation behavior  is  influenced  by  a  number  of  factors  such  as  the 
physical  structure,    porosity,    existence  and  movement  of  fluids  in  the 
pores,    etc.     In  addition,    such  geologic  features  as  faults,    joints,    seams, 
crushed  zones,    fissures  and  other  tectomic  effects  produce  behavior 
significantly  different  from  that  derived  on  the  assumption  of  continuous 
mass  [Ref.    14].     The  question  is  then  raised:     does  the  approximation 
based  on  the  Winkler  hypothesis  still  hold  over  various  types  of  soil 
foundations?     If  it  does  not,    then  how  much  error  may  be  expected? 
How  much  does  the  shearing  effect  in  the  material  of  the  foundation 
contribute  to  its  deformation?     Considering  only  the  deflections  of  one- 
dimensional  beams  on  foundations   (for  simplicity),    this  thesis  attempts 
to  answer  these  questions. 
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II.     DESCRIPTION  OF  ELASTIC  FOUNDATIONS 

This  chapter  presents  brief  descriptions  of  the  Winkler,    Modified 
Winkler  and  Continuous  foundations  considered  in  this  thesis.      Many- 
other  foundations  have  been  proposed  over  the  years  by  numerous 
investigators.     Kerr  [Ref.    31]  presents  a  summary  of  a  number  of 
linear  foundations. 

A.      DEFINITION  OF  WINKLER,    MODIFIED  WINKLER  AND 
CONTINOUS  FOUNDATIONS 

The  following  definitions  will  be  used  throughout  this  thesis. 

1.  Classical  Winkler  Foundation 

This  type  of  foundation  is  characterized  by  the  fact  that  the 
deflection  at  every  point  of  the  foundation  is  linearly  proportional  to  the 
pressure  applied  at  that  point  and  is  independent  of  pressures  or  deflec- 
tions acting  elsewhere  in  the  foundation.     This  assumption  is  equivalent 
to  considering  the  foundation  as  a  discontinuous  medium  composed  of 
independent  elastic  springs.     It  is  believed  that  this  hypothesis  was 
proposed  in  1867  by  Winkler     [Refs.    27,    9,    29]. 

2.  Modified  Winkler  Foundation 

This  type  of  foundation,    still  a  discontinuous  medium,    is 


Some  authors  claim  that  Euler  seems  to  have  been  the  first  to 
formulate  the  hypothesis,    although  it  is  generally  attributed  to  Winkler. 
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however  characterized  by  the  fact  that  the  pressure  at  every  point  is 
nonlinearly  proportional  to  the  deflection  at  that  point. 

In  this  thesis,    reference  to  a  Winkler  foundation  means  either 
the  classical  (linear)  Winkler  foundation  or  the  modified  (nonlinear) 
Winkler  foundation. 

3.  Continuous   Foundation 

In  contrast  to  the  Winkler  foundation,    the  pressure  and  deflec- 
tion at  a  point  in  a  continuous  foundation  is  affected  by  the  behavior  of 
the  entire  foundation  [Refs.    1,    9,    2?]. 

4.  Mixed  Mode  Foundation 

The  real  foundations  of  natural  soils  are  more  complex  than 
either  the  Winkler  or  continuous  foundations.      They  are  neither  com- 
pletely continuous  nor  discontinuous  [Refs.    14,    6].     The  pressure  and 
deflection  at  any  point  are  nonlinearly  dependent.     This  type  of  founda- 
tion is  therefore  considered  as  a  combination  of  the  modified  Winkler 
foundation  and  the  continuous  foundation. 

B.      MATHEMATICAL  MODELS 

In  this  chapter  the  governing  equations  of  the  foundations  considered 
in  this  thesis  are  developed. 

Consider  an  elastic  foundation  which  is  a  compressible  single  layer 
of  thickness  H  placed  on  a  rigid  base.     It  will  be  assumed  that  the  thick- 
ness of  the  elastic  foundation,    its  support  conditions,    the  elastic  con- 
straints and  all  other  properties  as  well  as  the  external  load  do  not 
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vary  in  the  z-direction.     A  narrow  plate  of  thickness  w  is  cut  from  the 
elastic  foundation  by  two  planes  parallel  to  the  x  y  plane.     A  beam  of 
the  same  thickness  w  rests  on  the  surface  of  this  elastic  foundation. 


lateral  loading 
q(x) 


Rigid  base 


Figure  1.      Geometry  of  an  elastically  supported  beam 


Let  an  external  load  q(x),  pounds  per  unit  length,  act  on  the  beam 
(Figure  1).  If  r(x)  is  the  reaction  of  the  elastic  foundation  against  the 
beam,    the  differential  equation  of  bending  for  the  beam  is  then  [Ref.    8]: 


(EIV")       =    q{x)  -   r  (x) 


(1) 
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where  EI  (x)  and  V  (x)  are  the  flexural  rigidity  and  the  deflection  of 
the  beam,    respectively. 

Equation  (1)  is  true  for  any  type  of  foundation.     It  contains  two 
unknown  functions,    V  (x)  and  r  (x).     In  order  to  determine  them,    an 
additional  relationship  is  needed.     This  relationship  is  associated  with 
the  response  of  the  elastic  foundation  and  depends  on  the  type  of 
foundation  being  considered. 

1.       Winkler  Foundations 

a.  Governing  Equation 

The  Winkler  foundation  is  characterized  by  a  foundation 
modulus,    generally  denoted  by  the  letter  kt     Let  V,    inches  and  r,    pounds 
per  inch,    be  the  deflection  and  the  reaction  of  the  foundation,    respectively. 
The  Winkler  foundation  is  characterized  by  the  equation 

r  (x)  =  k*vP(x)  (2) 

For  the  classical  Winkler  foundations,     p  =  1  and  k*is  a 
constant  and  is  expressed  in  pounds  per  cubic  inch.      This  assumption 
has  proved  adequate  for  calculating  stresses  and  deflection  in  railroad 
tracks  [Ref.    1]. 

b.  Values  of  Foundation  Modulii 

For  sand,    k*may  vary  from  25  to  100.    Ib/cu.  in.      For 
ordinary  soils  on  which  railroad  tracks  have  been  built,    k  may  vary 
from  110.    to  130.    Ib/cu.  in.      For  gravel,    the  value  of  l^is  even  more 
uncertain  and  it  can  vary  from  200.    to  1200.    Ib/cu.  in.    [Refs.    1,    27]. 


17 


2,  Modified  Winkler  Foundations 

For  actual  foundations  of  natural  soils,    k*is  the  coefficient 
associated  with  the  nonlinear  term  [Ref.    6].      Therefore,    to  generalize 
the  Winkler  hypothesis,    let  p    be  any  positive  real  number.     Such 
foundations  are  called  modified  Winkler  foundations. 

3.  Continuous   Foundations 

a.  Governing  Equation 

Let  s  be  the  shear  force  in  the  foundation,    pounds  per 
inch,    and  l^be  the  foundation  modulus,    pounds  per  cubic  inch.      The 
shear  force  s  is  a  measure  of  the  friction  force  between  soil  particles. 
The  governing  equation  of  continuous  foundations  is  [Ref.    9]: 

r  (x)  =  k*V  (x)  -  (s  v'(x))'  (3) 

Details  of  the  continuous  foundation  analysis  are  given  in 
Appendix  A. 

b.  Values  of  s  and  k 

The  values  of  s  and  k  depend  on  the  contact  area  of  beam 
and  foundation  per  unit  length  of  beam.      For  a  loaded  area  of  6.  5  inches 
of  width,    equal  to  that  of  a   12WF27  standard  beam,    and  1.    inch  of  length, 
typical  values  of  s  and  k  for  some  soils  are: 

Average  sea  floor  sediment: 

s  =  18.6  lb  k  =  8.  3  Ib/in^ 

A  typical  beach  soil: 

s  =  1000.    lb  k  =  500.    Ib/in^ 


where  k  =  wk ' 
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A  typical  inland  soil: 

s  =  6000.    lb  k  =  3000.    lb/in 

Additional  values  of  soil  characteristics  are  given  in 

Reference  7. 

4.       Mixed  Mode  Foundations 
a.       Governing  Equation 

The  governing  equation  of  this  type  of  foundations  can  be 

derived  from  the  definition  of  the  mixed  mode  foundations  given  in 

part  A.      The  mixed  mode  foundation  is  a  combination  of  a  modified 

Winkler  foundation  and  a  continuous  foundation: 

r  (x)  =  k  vP(x)  -  (s  V  (x))'  (4) 

It  is  seen  that  any  of  the  previous  foundations  are  special 

cases  of  the  mixed  mode  foundation.      For  the  Winkler  foundation  s  =  0 

and    p   =   1,    for  the  nonlinear  foundation  s  =  0  and  for  the  continuous 

foundation    P   =  1. 
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III.  GENERAL  FORMULATION  OF  THE  PROBLEM 

A._  GOVERNING  EQUATIONS 

Consider  a  beam  on  an  elastic  foundation.     The  differential 
equation  of  bending  for  a  beam  on  an  elastic  foundation  is  given  by 
[Ref.    8]: 

(Eiv"(x))"     =    q(x)  -   r(x)  (1) 

where  q(x)  is  the  external  loading,    V(x)  is  the  beam  deflection  and  r(x) 
is  the  reaction  of  the  foundation.      The  mixed  mode  foundation,    by- 
definition,    can  be  considered  as  the  most  general  type  of  foundation. 
Its  governing  equation  is  given  Eqn.    1  with 

r(x)  =  kvP  (x)  -  (sv'(x))'  (4) 

and  p  >  0. 

Substitution  of  (4)  into  (1)  gives: 

(EIV    (x))"     =    q(x)  -  kvP   (x)  +  (sv'(x))  (5) 

or 

(EIV)      _   (sV)     +  kV^'    =  q  (6) 

To  solve  this   equation  for  any  positive  real    p   ,    consider  a  least 
square  best  fit  polynomial  such  that 

V^(x)  =  bj  +  b^  V(x)  +  b^  V^(x)  +  .  .  .  (7) 

where  the  b.  coefficients  (i  =   1,2,3,...)  are  constants. 


20 


For  simplicity  in  the  development  of  an  actual  computer  program, 
equation  (7)  is  taken  to  be  a  second  degree  polynomial.     This  restriction 
is  not  inherent  in  the  Galerkin  procedure  but  rather  in  the  Finite  Element 
formulation  of  Galerkin  adopted  in  this  thesis.     Substitution  of  (7)  into 
(6)  leads  to  the  following  equation: 

(EIV)"  -  (sV)'  +  k  (b^  +  b2V  +  b^V^)  =  q  (8) 

or 

(EIV")"  -  (sV)'  +  kb^V  +  kb^V^  =  q  -  kbj  .  (9) 

Letting 

a        =  kb^  (10) 

^2      =  kb^  (11) 

f  =  q(x)   -   kbi  (12) 

equation  (9)  may  be  rewritten: 

(EIV")"  -   (sV)'  +  a^V    +  a^V       =    f  (13) 

1.       Beam  on  a  Winkler   Foundation 


Equation  (13)  can  be  applied  to  the  case  of  a  beam  on  a  classical 
Winkler  foundation,    by  letting  bj^  =  b^  =  0,   b2  =  1  and  Sj,    =  0.     Hence 
the  governing  equation  of  bending  for  a  beam  on  a  classical  Winkler 
foundation  is: 

(EIV")"  +  ajV  =  q(x)  (14) 

where  a,   is  defined  by  Eqn.    10. 

2.       Beam  on  a  Modified  Winkler  F'oundation 

When  s  =  0,    Equation  (13)  becomes  the  governing  equation 
of  bending  for  a  beam  on  a  modified  Winkler  foundation: 
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(EIV")"  +  a^V  +  a^V^     =    f  (15) 

where  terms  are  defined  by  Eqns.    10,    11,    12. 

3.  Beam  on  a  Continuous   Foundation 

Equation  (13)  can  be  applied  to  the  case  of  bending  of  a  beam 
on  a  continuous  foundation,    by  letting  b     =  1,    b,    =  b-   =  0  : 

(EIV")"  -  (sV)'  +  a^V  =  q(x)  (16) 

4.  Beam  on  a  Mixed  Mode  Foundation 

Equation  (13)  is  the  governing  equation  of  bending  for  a  beam 
on  a  mixed  mode  foundation: 

(EIV")"  -  (sV)'  +  a^V  +  q^V^  =  i  (17) 

B.  BOUNDARY  CONDITIONS 

Reference  8  thoroughly  develops  the  boundary  conditions  for  a  beam. 
Only  the  boundary  conditions  on  displacement  and  slope  (the  so  called 
essential  or  principal  boundary  conditions)  must  be  imposed  in  a 
variational  formulation  of  the  problem. 

C.  ASSUMPTIONS  AND  RESTRICTIONS 

Some  restrictions  have  already  been  stated  in  the  introduction 
section  of  this  thesis;    other  restrictions  will  now  be  discussed. 

1.       The  flexural  rigidity  EI  must  be  "slowly  varying"  in  accord- 
ance with  the  development  of  the  basic  equations  of  classical  beam 
theory. 
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2.  The  shear  coefficient  s  must  also  be  "slowly  varying"  in 
accordance  with  the  development  in  Appendix  A. 

3.  The  foundation  modulus  k(x)  and  the  load  variable  q(x)  need  not 
be  slowly  varying. 

D.      SOLUTION  TECHNIQUES 

The  Finite  Element  Method  via  the  Galerkin  approach  will  be  used 
here.     Its  formulation  will  be  presented  in  the  next  section.      There  are 
two  types  of  problems  to  be  considered. 

1.  Linear  problems 
Consider  the  linear  system 

(EIV")"  -  (sV)'  +  kV  =  q  (18) 

Suppose  the  solution  for  q  =  q     is  V,   and  the  solution  for  q  = 

^2  ^^  ^7'    then  for  q  =  qi   +     °-  ^^  we  have  the  solution: 

V  =  Vi   +    a  V^  (19) 

The  application  of  the  Finite  Element  Method  to  any  linear 

differential  equation  leads  to  a  linear  system  of  algebraic  equations 

which  will  be  solved  by  Gaussian  elimination. 

2.  Nonlinear  problems 

In  contrast  to  the  linear  case,    any  nonlinear  differential 
equation,    like  (13)  or  (15),    leads  to  a  nonlinear  system  of  algebraic 
equations  which,    in  this  thesis,    will  be  solved  by  Brown's  iteration 
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method.     Reference   13  thoroughly  develops  the  Brown's  iteration  method 
which  is  a  modified  Newton-Raphson  method. 

The  results  of  the  Finite  Element  Method  solutions  presented 
here  will  be  checked  with  the  classical  solutions  if  available,    or  with 
other  methods,    such  as  the  Finite  Difference  Method.      Programming 
details  will  be  presented  in  Chapter  V. 
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IV.      FINITE  ELEMENT  FORMULATION 

A.     APPROXIMATION  BASED  ON  THE  GALERKIN  PROCEDURE 

In  order  to  obtain  a  nixmericaJ  solution  for  displacements  of  beams 
on  elastic  foundations,   let  V  be  approximated  by  the  m  degrees  of 
freedom  expression: 

m 

V    «      XI    ^i  "^i  (^^  (2^) 

i=  1 

where  G.(x)  are  global  (or  system)  shape  functions  and  the  unknowns 

V-  coefficients  are  generalized  coordinates  [Ref.    14],    i.e.,    system 

degrees  of  freedom.     Since  the  Galerkin  method  deals  directly  with  the 

differential  equation  [Refs.    21,    15,    22],    (21)  will  be  inserted  in  (13) 

and  the  equation  residual  is  formed  as  follows: 

m  ,,  , 

J2      [<(EIV.G.    )",    G^(x)  >  -      <   (s        V.G.'(x)).    Gj^(x)> 

i=  1 

m^ 
+    <ai     V.G.(x),    Gi^(x)  >       +   XI    <S       V.   G.   (x)y.G  (x),G^(x)    >    ] 

i=i  ^     ^         ^   ^ 

-     <   f.    Gj^(x)    >       =0 
where  the  notation  <  a,b>        means  lQa(x)b(x)dx     .    Then 

J2    [V.      <    (EIG!'(x))",    Gj^(x)    >       -  V.     <  (s  g!(x)),    Gj^(x)    > 
1=1 

m 
+  a^V.    <     G.(x).    Gj^(x)  >      +    J2      ^Z^^j    <     Gi(x)Gj(x).  Gj^(x)   >      ] 

-     <  f,  G.  (x)    >       =  0 
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Integration  of  the  first  two  terms  by  parts  yields: 


G^  (X)    > 


J]     [V.     <    EIG."(x).    GUx)  >  +      V.     <     sG.'(x), 

1=1 

m 
+  a^V.    <    G.(x),    Gj^(x)>      +   ^       ^z'^i'^j    ^    G.(x)G.(x),G^(x)    >      ] 

-     <    f,    Gj^(x)>      =0  (22) 

There  is  a  constant  of  integration  due  to  boundary  conditions  left 
by  the  process  of  integration  by  parts  in  Eqn.    22.     However  this  constant 
can  be  omitted  in  the  Galerkin  procedure  [Ref.    23]. 

Let 

K..  =  <EIg!'(x),    GJ^(x)> 

<;sg!(x),    G^(x)> 

ik  =  <  Gi(-)'    ^k^-)  > 

ijk  =  <G.{x)G.(x).    Gj^(x)> 

Qk  =  <f.Gj^(x)> 


ik 
Cik 


M 


N, 


(23) 


Then,    Eqn.    22  may  be  rewritten  as  follows: 


m 
C        [K.^  +  s  Cik+a^Mik]     V. 


i=l 


m  m 


i=i   j=i 


(24) 


for  k=l,2,3,...,m 

This  is  a  system  of  nonlinear  algebraic  equations,    where  the 
unknown  displacements  V.  must  also  satisfy  the  imposed  boundary 
conditions. 
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B.      ELEMENT  SPiAPE  FUNCTIONS 

Consider  a  beam  element.     There  are  two  nodal  coordinates,    the 
deflection  and  slope  at  each  end,    and  therefore  each  element  must  have 
a- total  of  four  degrees  of  freedom.     On  the  element  level  these  coor- 
dinates are  numbered  as  follows: 

Coordinate  1      =     V(o)  Coordinate  3      =     V(i) 

t  t 

Coordinate  2      =     V  (o)  Coordinate  4      =     V  (t) 

The  compatible  element  shape  functions  must  be  cubics  and  it  is 

evident  that  four  independent  shape  functions  are  required  [Refs.    14,    16]. 

e 

Let  G.  (x),    the  element  shape  functions,  be  defined  as  follows,    where  the 

superscript  "e"  refers  to  the  beam  element: 
i 


V(o) 


g:  =  i 


ci     2  x      7      3 

3x  +    Zx 


lr3 


Cj~    =        X  _     Z      X         +      X 


(25) 


V(i)       G^     =      Sx'^       -     2 


G' 


l2 

3 


Figure  2.     Element  shape  functions 


Note  that  the  function  G.(x)   represents  the  variation  of  V  along  the 
element  in  such  a  way  that  g|(o)  =  G|'(o)  =  G®(1)  =  G|'(i)  =  1  and  that 
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other  evaluations  of  G   (o),    G.  (1),    G       (o),    G.     (1)  vanish;  cf.    the  sketches 

i  1  1  1 

in  figure  2.    where  x  is  the  local  coordinate  from  left  to  right  and  1 
denotes  the  length  of  the  beam  element. 

C.      FORMATION  OF  ELEMENT  MATRICES 

Once  the  element  shape  functions  are  defined,    the  components  of 

e  €         e  e  g 

the  element  matrices    K.,  ,    C,  ,    M.,  ,    Q,      and  N  can  be  computed 

ik       ik         ik        ^  ijk 

through  ordinary  integration  of  polynomials.      The  following  results, 
associated  with  the  linear  terms,   have  been  given  in  many  text  books 
[Refs.    21,    14,    25]. 

12.  6.       -      12.  6. 


K^^    =     <EIG^^    (x),    G^   (x)>  =    El' 


i3             i2 

i3 

i2 

4. 
1 

-6. 
l2 

2. 
4 

Symmetric 

12. 
t3 

-6. 
1^ 

4. 
ir 

where  EI     is  a  constant  (the  average  EI)  over  each  element. 


e  >     e'        e'       . 

C        =<;sG.Cx),Gj^  (x)'>=     s* 


ik 


1.2  0.  1  -1.2 

1  I 

.  1333331    -0.  1 


symmetric 


1.  2 


0.1 

-0.0333331 
_0.1 
0.133333  1 


where  s     is  a  constant  (the  average)  over  each  element. 
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<k  = 


<G^^(x),  G^x)>- 


Q^=    <f.G^(x)>=      f 


2  2 

.3714291:       .0523811        .1285715:     -.03095^1- 

.0095241      .0309521r        -.007143]r^ 


symmetric 


37l4291r       -.0523811 


.0095241 


.51 

2 

. 0833331 

.  51 
-.083333l' 


The  components  of  N. .,     associated  with  the  nonlinear  response  are 
not  in  the  literature.     They  are  as  follows: 


N 
N 
N 
N 
N 
N 


e 

ijk 
e 

111 
e 

211 

e 

311 

e 

411 

e 
221 


321 


N 
N 
N 


421 

e 

331 

e 

431 


<G.   (x)  G^(x),    G^(x)> 

. 3071431 

2 
.0384921 

.0642861 

2 

-.0170601 

. 0063491^ 
.0138891^ 


0035711- 


.0642861 


-.0138891 
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<! 

.0031751 

^222 

.0011911'^ 

^322 

.0031751 

^^22 

4 
-.0007941 

^3%2 

2 

.0170641 

^132          ^ 

-.0035721^ 

^42          ^ 

4 
.0007941 

^333          ^ 

.3071431 

^^433 

-.0384921^ 

^^43          ^ 

.0063491^ 

^^44 

4 

-.0011901 

e 
Notice  that,    according  to  Eqn.    23,    the  N. .,    coefficients  are  equal 

e  e  e 

for  any  permutation  of  i,    j,    k  (i.  e.  ,    N112  ~  ^121  ~  •'^211'    ^^^'  ) 


D.      FORMATION  OF  SYSTEM  MATRICES 

The  system  matrices  are  obtained  from  the  element  matrices  as 
follows : 

1.  A  correspondence  table  between  local  and  global  coordinates 
is  formed 

2.  The  element  coefficients  of  any  matrix,    say  the  element 

e 

stiffness  coefficients  K.  •,    are  substituted  directly  into  the 

system  stiffness  matrix  K_    according  to  the  correspondence 
of  element  coordinates  i,    j  to  system  coordinates  I,    J.[Ref.l4j 
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V.     COMPUTER  PROGRAMS 

Equation  (13)  is  originally  of  the  form: 

(EIV")"  -  (SV)'  +  kvP         :r  q  -  (29) 

The  following  block  diagram  shows  how  the  investigation  proceeds 
(Figure  3).      Equation  13  or  29  will  be  solved  by  the  Finite  Element 
Method.      First  consider  the  case  of  a  classical  Winkler  foundation, 
for  which  the  theoretical  analysis  is  readily  available.      The  results 
from  the  Finite  Element  Method  will  then  be  compared  to  the  theoretical 
results.     Similarly,    the  modified  Winkler  foundation  will  be  solved,   also 
by  the  Finite  Element  Method.     Since  there  is  no  theoretical  analysis 
for  the  case  of  these  nonlinear  problems,    some  of  the  results  from  the 
Finite  Element  Method  will  be  compared  to  those  from  the  Finite  Dif- 
ference Method.      Finally,    in  dealing  with  continuous  foundations,    cer- 
tain types  of  soils  will  be  considered.      Results  obtained  by  varying 
flexural  rigidity  of  beams,    foundation  modulus  and  applied  force  will  be 
presented  in  the  next  chapter. 

A.   GENERAL  PROGRAM  ORGANIZATION 

The  program  may  be  divided  into  two  separate  categories,    depend- 
ing on  linearity  or  nonlinearity.      The  reason  for  this  distinction  is  that 
for  nonlinear  problems,    the  Brown's  iteration  technique  is  applied  and 
double  precision  is  required.     In  contrast,    a  simple  Gaussian  elimina- 
tion method  may  be  applied  to  solve  the  linear  systems. 
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(EI  V")"  .   (s   V  )'   +  kvP  =  q 

'' 

t 

\i 

Winkler  Foundation 
p  =  1.      s  =  0 

Modified 
Winkler 
Foundation 

Mixed  Mode 
Foundation 

• 

^ 

\f 

> 

f 

\f 

\' 

Theoretical 
Solution 

Finite 

Difference 

Method 

Finite  Element 
Method 

\ 

1 

f 

Results  Compared 

\ 

t 

\ 

/ 

Four 
Elements 

Eight 
Elements 

Variable  Parameters 
EI ,    k  ,    q  ,    s 

/ 

\                   - 

> 

f 

> 

i 

i 

I 

Results 

Tabulj 

It 

ed 

Test  for 
Convergence 

Figure  3.     Sequence  of  the  problem  investigation 
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( 


START 


) 


(1) 

Store  all  informa- 
tion needed  for  the 
problem.      (Number  of 
elements  is  also 
specified) 


-> 


JiL 


(2) 


Supply 
necessary  data 


^ 


(3)  Matrix  formation 


J±. 


(5) 


Actual  problem 
solution 


NO 


■> 


(4)   Finite  Diff.    Method 
solution  (optional) 

See  Figure  7 


YES 


< 


\^ 


STOP 


) 


Figure  4.     General. organization  of  the  program 
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1 — 

'     (3) 

STIFFl 

/       STIFF2 

\^ 

LOAD 

N 

/ 

BOUND 

Figure  5.     Matrix  formation  (STIFFl  is  used  for  classical 
Winkler  foundation  only,    STIFF2  is  used  for  all  other  cases.  ) 
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r 


(4) 


SOLVE 


\L 


RESULT 


ZSYSTM 


Jk 


RESULT 


4- 


AUX 


(a) 


(b) 


J 
J 

J 


Figure  6.     Actual  problem  solution: 

(a)  for  linear  problem, 

(b)  for  nonlinear  problem. 
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(4) 


1 


DIFSOL 


±. 


ZSYSTM 


±. 


RESULT 


■<r 


DIFF4     /     DIFF8 


L 


Figure  7.      Finite  Difference  Method  solution 
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Each  type  of  program  performs  four  major  functions: 

1.  Storing  all  information  needed  for  subsequent  problems 

2,  Supplying  necessary  data  to  each  particular  problem 
3. a.    Forming  the  element  matrices 

b.  Forming  the  system  matrices 
4.       Solving  the  algebraic  problem. 
This  organization  is  shown  in  Figures  4  through  6. 

The  Finite  Difference  Method  is  applied  to  some  problems  to  check 
the  results.     Its  general  flow  chart  is  shown  in  Figure  7. 

B.     SUBROUTINES 

The  MAIN  program  is  a  control  segment  that  serves  to  call  other 
subroutines  in  proper  order  as  required.     No  calculations  are  performed 
in  the  MAIN  program  (See  Appendix  D). 

1.  Subroutine  STORE 

This  subroutine  reads  all  data  available  for  the  problem  under 
investigation,    including  the  maximum  number  of  elements,    the  number 
of  nodes,    the  number  of  boundary  conditions,    the  number  of  degrees  of 
freedom,    the  correspondence  table  of  global  and  local  nodes,   load 
conditions,   beam  properties  and  foundation  characteristics. 

2.  Subroutine  TEST 

This  subroutine  is  used  particularly  for  testing  the  convergence 
of  the  solution.     In  doing  so,    it  selects  special  information  and  supplies 
it  to  the  MAIN  program  which,    keeping  all  information  constant,    doubles 
the  number  of  elements  from  one  run  to  another. 
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3.  Subroutine  CURFIT 

For  the  case  of  a  nonlinear  elastic  foundation  for  which  the 
reaction  r(x)  =  kV         ,    where    p    is  not  an  integer,    a  direct  application 
o"f  the  finite  element  method  via  the  Galerkin  procedure  is  not  efficient. 
Using  a  least  square  best  fit,    this  subroutine  replaces  kV^      by  a 
second  order  polynomial. 

4.  Subroutine  INCHK 

Depending  on  various  conditions  in  changing  beam  flexural 
rigidity  and  foundation  modulus,    this  subroutine  selects  information 
from  subroutine  STORE  and  supplies  it  to  the  MAIN  program.      This 
subroutine  also  prints  the  appropriate  echo  check. 

5.  Subroutines  STIFFl  and  STIFF2 

Subroutine  STIFFl  is  used  for  the  linear  problem  and  uses 
single  precision;  subroutine  STIFF2  is  used  for  the  nonlinear  problem 
and  uses  double  precision.      Both  of  these  subroutines  form  the  element 
stiffness  matrices  and  assembles  the  system  stiffness  matrix. 

6.  Subroutine  LOAD 

This  subroutine  forms  the  element  load  vector  and  assembles 
the  system  load  vector.     It  can  be  used  for  both  linear  or  nonlinear 
problems,    provided  it  is  in  double  precision  for  the  later  case. 

7.  Subroutines  BOUNDl  and  BOUND2 

These  subroutines  apply  the  boundary  conditions  to  the  system 
force  vector  and  the  system  stiffness  matrix  produced  by  subroutines 
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LOAD,    STIFFl  or  STIFF2,    and  transform  them  into  an  appropriate 

form  ready  to  be  solved.     Subroutine  EOUNDl  is  used  for  linear  problems 

and  subroutine  BOUND2  for  nonlinear  problems. 

8.  Subroutine  SOLVE 

This  subroutine  employs  the  Gaussian  elimination  method  to 
solve  the  linear  system  of  algebraic  equations  for  the  linear  problem. 

9.  Subroutine  RESULT 

This  subroutine  prints  any  output  data  from  either  linear  or 
nonlinear  problems.      For  nonlinear  problems,    double  precision  must 
be  Supecified. 

10.  Subroutine  DIFSOL 

This  subroutine  solves  the  nonlinear  problem  by  the  Finite 
Difference  Method  of  which  the  development  will  be  presented  in 
Chapter  V,    part  D.      This   solution  is  optional  and  was  employed  in  a 
few  cases  in  order  to  check  the  results  of  the  nonlinear  problem  given 
by  the  Finite  Element  Method. 

11.  Subroutine  FSOIL 

Similar  to  subroutine  CURFIT,   using  a  least  square  fit,    this 
subroutine  replaces  the  total  settlement  curve  of  a  natural  soil  given 
by  experimental  data  [Ref.    30],    by  a  polynomial  of  second  degree.      The 
total  settlement  curve  of  the  sea  floor  sediment  is  supposed  to  have  the 
same  shape  as  that  of  the  inland  soil  but  its  magnitude  is  reduced  by  a 
certain  scale. 
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12.  Subroutine  VGUESS 

This  subroutine  is  designed  to  give  an  initial  estimate  vector 
to  the  nonlinear  algebraic  system  before  any  iteration  process  is  made. 
It  is  desirable  that  if  the  first  estimate  has  failed  to  make  the  iteration 
convergent,    another  estimate  is  provided.     In  this  subroutine,    the 
initial  estimate  must  be  taken  in  the  form  of  a  second  order  polynomial. 

13.  Function  AUX 

This  is  a  nonlinear  algebraic  equation  system  translated  from 
equation  (24),    which  has  to  be  solved  and  which  serves  as  an  external 
input  to  subroutine  ZSYSTM. 

14.  Functions  DIFF4  and  DIFF8 

These  are  nonlinear  algebraic  equation  systems  translated 
from  Eqns.  27  and  28,  respectively,  which  have  to  be  solved  by  the 
Finite  Difference  Method  and  which  serve  as  external  inputs  to  sub- 
routine DIFSOL. 

15.  Subroutine  LSQPL2 

This  library  subroutine,    from  the  Naval  Postgraduate  School, 
Monterey,    California,    employs  a  least  square  best  fit  method  to  replace 
any  function  by  a  polynomial. 

16.  Subroutine  ZSYSTM 

This  library  subroutine,    from  IMSL,   solves  nonlinear 
algebraic  equations  employing  the  Brown's  method  of  iteration,    which 
is  a  modified  Newton's  method  and  is  developed  in  Reference  13. 
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C.      TESTS  FOR  CONVERGENCE 

The  test  problems  employed  here  are  designed  to  verify  the  opera- 
tion of  the  computer  program  and  test  the  convergence  of  the  equation 
system.     The  two  categories  of  test  performed  were  linear  and  non- 
linear problems. 

An  analytic  solution  and  the  Finite  Difference  Method  solution  were 
used  to  verify  the  Finite  Element  Method  solution  for  linear  and  non- 
linear cases  respectively.     The  following  results  (Tables   1,    2,    3)  show 
that  for  an  uniform  beam  of  100  inches  in  length,    with  a  uniform  load, 
the  8-element  model  gives  sufficiently  accurate  results.  These  analyses 
show  that  as  the  number  of  elements  is  increased  the  finite  element 
method  results  approach  the  theoretical  results.      The  small  difference 
between  eight  and  sixteen  element  results  (less  than  1%)  suggests  the 
solution  has  converged  sufficiently  for  engineering  purposes.     For 
symmetric  problems  (i.  e.  ,    symmetric  EI,    k,    s  and  q),    advantage 
may  be  taken  of  syrametry.     In  this  case,    only  half  the  system  need  be 
considered,    thus  greatly  reducing  computer  effort.      This  reduction 
is  especially  beneficial  for  nonlinear  problems.    The  time  consumption 
may  increase  from  10  to  30  fold  in  going  from  a  4-element  model  to  an 
8-element  model,    depending  on  how  good  the  initial  guess  is.     A  brief 
description  of  the  test  program  is  shown  in  Figure  8. 
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N^ 


) 


Problem  counter 


Store  all  information 
needed  for  the  prob- 
lem.     (Number  of 
elements  is  also 
specified) 
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■> 
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I 
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of  elements 
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STOP 


) 


Figure  8.      Flow  chart  of  a  test  program 
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)  EI  =  S.xlO*^  Ibs-in 

(  q     =  22.  5  lbs/in 


Figure  9.      Compared  deflection  curves  of  a 
typical  inland  soil  from  Figure  15, 
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Deflection 
in  inches 


0.  8 


O     Sea  floor  model  foundation 
A     Winkler  foundation 


1.       2.25     4.5      6.75      9 


15.75  22.5 

Load  lbs /in 


Figure  10.      Deflection  of  a  beam  on  a  sea  floor  model  foundation 
compared  to  that  of  Winkler  foundation. 
(Maximum  difference  is  about  6%) 
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Table  1.      Convergence  of  the  Linear  System 
(EIV")"   +  kV  =  q 
Uniform  Simply  Supported  Beam  with  Uniform  Load 
L  =   100  in,    EI  =  3  X  10^  Ib-in^, 
q  =  1  lb/in,    k  =   1  Ib/in^, 
Deflection  Given  in  Inches  ( /3L  =  0.  96) 


Location 

X 

/L 

0. 

0 

0. 

125 

0. 

250 

0. 

375 

0. 

,500 

Finite  Element  Method  Solution 
Number  of  Elements        Number  of  Elements 
8  16 


0.0 

0.016296 

0.029897 

0.038839 

0.41950 


0.0 

0.016310 

0.029919 

0.038863 

0.041975 


Theoretical 
Solution 


0.0 

0. 016301 

O.O299O6 

0.038880 

0.042178 


Table  2.     Convergence  of  the  Nonlinear  System 
(eTv")"  +  kV^   =:  q 
Uniform  Simply  Supported  Beam  with  Uniform  Load 
L  :=   100  in,    EI  =  3  X  10''  Ib-in^ 

q  =  0.  01  lb/in,    k  =   1  Ib/in^ 
Deflection  Given  in  Inches  (   /3L  =  0.  96) 

Finite  Element  Method  Solution: 


Location 
x/L 


Number  of  Elements 
4 


Number  of  Elements 
8 


0.0 

0.250 

0.500 


0.0 

0.  1451D-03 

0. 1483D-03 


0.0 

0.0456D-03 
0.  1485D-03 
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..  2 

Table  3.     Solution  of  (EIV")     +  kV     ^  q 

F,  E.  M.    Compared  to   F.  D.  M.    Uniform 

Simply  Supported  Beam  with  Uniforrn  Load 

L  =  100  in,    EI  =  6.  123x10^  Ib-in^, 

q  =  2.25  lb/in,    Number  of  Element  =  8 


Location 

k  =  100 

/3L  =  .80 

x/L 

F.  E.M. 

F.  D.M. 

ITMAX  =  3 

ITKIAX  =  3 

0.  0 

0.0 

0.0 

0.  125 

0.  185781D-03 

0.  188397D-03 

0.250 

0. 340909D-03 

0.  345395D-03 

0.  375 

0.442958D-03 

0.448565D-03 

0.500 

0.478469D-03 

0.  484450D-03 

Location 

k  =  500 

^L  =  1.20 

x/L 

F.  E.M. 

F.  D.  M. 

ITMAX  =  8 

ITMAX  =  3 

0.0 

0.0 

0.0 

0.  125 

0.  185776D-03 

0.  188392D-03 

0.250 

0.340900D-03 

0.  345385D-03 

0.375 

0.442946D-03 

0. 448552D-03 

0.500 

0.478456D-03 

0. 484436D-03 

jocation 

k  =  1000 

^L  =  1.42 

x/L 

F.  E.M. 

F.  D.M. 

ITMAX  =  4 

ITMAX  =  3 

0.0 

0.0 

0.0" 

0.  125 

0.  185769D-03 

0.  188385D-03 

0.250 

0.  340889D-03 

0.  345373D-03 

0.  375 

0.442931D-03 

0.448536D_03 

0.500 

0.478440D-03 

0. 484419D-03 
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Location 

k  =  2000 

j3L  =  1.69 

x/L 

F.  E.M. 

F.  D.M. 

ITMAX  =  4 

ITMAX  =  3 

0.0 

0.0 

0.0 

0.  125 

0.  185757D-03 

0.  188372D-03 

0.250 

0.340865D-03 

0.  345349D-03 

0.375 

0.442900D-03 

0.448504D-03 

0.500 

0.478407D-03 

0.484384D-03 

Location 

k  =  3000 

^L  =  0.  87 

x/L 

F.E.M. 

F.  D.M. 

ITMAX  =  4 

ITMAX  =  4 

0.0 

0.0 

0.0 

0.  125 

0.  185745D. 

.03 

0.  188359D-03 

0.250 

0.340842D- 

03 

0.345324D-03 

0.  375 

0.442870D- 

03 

0. 448472D-03 

0.500 

0.478374D- 

03 

0. 484350D-03 
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D.      THE  FINITE  DIFFERENCE  METHOD 

Since  there  is  no  analytic  solution  to  the  nonlinear  case,    the  Finite 

Difference  Method  was  used  in  this  thesis  to  verify  the  results  of  the 

Finite  Element  Method.      The  development  is  restricted  to  uniform  EI 

and  k.      The  finite  difference  approximations  are  developed  in  Ref.    15. 

The  fourth  order  derivative  has  the  central  difference  computational 

molecule  as  follows: 

(±JL).  =J_  [V(j-2)  _  4V(j-l)+6V(j)-4V(j  +  l)+V(j+2)]  +0(1^) 
dx4     J         l4 

where  1  is  an  elemient  length 

Consider  the  differential  equation  of  an  uniform  simply  supported 

beam  with  uniform  load  on  a  uniform  nonlinear  foundation: 

EIV"  +  kV^  =  q 

or  EIv"  +    ^^   (EIV)     =  q 

-  k 

Let  V     =  EIV  and  ^     =    ^       >   then 

the  above  equation  may  be  rewritten  as  follows: 

V*"  +  A    V''^^  -  q  =  0 

A  computational  molecule  model  applied  at  node  j  yields: 

J_         v'''(j-2)  .  4V'^j-l)  +  6V'''(J)  _  4V*(j  +  l)+V*(j+2) 
X4 

+  A    Y'^'^ii)  -  q(j)  =  0  -         - 


or: 


V*(j-2)_4V""(J-1)  +  6V''(J)  :  4V'"'(J  +  1)  +  V*(j+2) 

A    ^'<2  4  <26) 

+  A    X^V      (j)  -  X   q(j)  =  0 
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where  j  ranges  from  1  to  2  for  4-element  model  and  from  1  to  4  for 

8- element  model. 

A.      FOUR- ELEMENT  MODEL 


V(-l)  =  -V(l)         V{0) 


V(l) 


V(2) 


V(l) 


• 


V(0) 

• 


X 


^v^d 


6  V(l)  -  4  V(2)  +  A   X   V^{1)  -  qX     =  0 
-8  V(l)  +  6  V(2)  +  A  x'^V^(2)  _  qX^  =  0 


(27] 


B.      EIGHT-ELEMENT  MODEL 

-V(l)        V(0)        V(l)        V(2)        V(3)        V(4)        V{3)        V(2)        V(l)        V(0) 

• • • m  • • • • • • 


{      (28) 


X 
5V(1)  -  4  V(2)  +  V(3)  +  A   X   V^(l)  _  qx"^  =  0 
-4V(1)  +  6V(2)  _  4V(3)  +  V(4)  +  A   x'^V^(2)  _  qX     =  0 
V{1)  -  4V(2)  +  7V(3)  -  4V(4)  +  A   X^V^(3)  -  qX'^  =  0 
2V(2)  _  8V(3)  +  6V(4)  +  A   X^V^(4)  -  qX     =0 


Notice  that  the  superscript  ■■'•  has  been  omitted  for  simplicity;    the 
result  V(j)  must  be  multiplied  by  (EI)"     to  get  the  actual  deflection. 

The  above  development  is  restricted  to  a  uniform  beam  on  a  founda- 
tion with  uniform  k.     In  this  case  the  finite  difference  method  has 
distinct  advantages  over  the  finite  element  method.     In  the  case  of 
variable  EI,    k,    and  s  the  finite  difference  method  is  not  easily  developed 
and  the  advantage  shifts  to  the  finite  element  method. 
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VI.     ILLUSTRATIVE  PROBLEMS  AND  RESULTS 

This  chapter  presents  the  results  of  two  investigations. 

The  first  investigation  attempts  to  compare  the  behavior  of  three 
foundation  models,    the  Winkler,    the  modified  Winkler  and  the  mixed 
mode  foundations.     The  comparison  is  limited  to  the  case  of  a  simply 
supported  uniform  beam  on  foundations  with  uniform  k  and  s  properties. 
Table  4  gives  values  for  the  parameters  considered  in  this  study.      The 
value  of  ^L  varied  from  1.  62  to  3.  98  i.  e.  ,   beams  of  medium  length. 
The  results  obtained  apply  only  to  the  particular  simply  supported 
systems  considered.     Other  problems  may  yield  other  results. 

The  second  investigation  seeks  to  establish  that  the  Finite  Element 
Program  developed  in  this  thesis  can  solve  problems  with  variable  EI, 
k,    s  and  q.     As  a  sample  problem,    the  case  of  a  free-free  beam  with 
variable  stiffness,    variable  foundation  modulus,    and  variable  load  was 
considered.      The  j3L  for  this  illustrative  problem  was  3.  5. 

A.      UNIFORM  SIMPLY  SUPPORTED  BEAMS  WITH  UNIFORM  LOADS 
ON  SOIL  FOUNDATIONS 

Figures  9  and   10  present  the  deflection  of  a  simply  supported  beam 

on  soil  foundations  which  are  either  inland  soil  or  sea  floor  sediment. 

The  range  of  parameters  considered  is  shown  in  Table  4  below. 
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Table  4.     Information  Data 


p 

Loading 
Conditions 

q 

(lb/in) 

Beam 
Properties 

^^2 
(lb- in'') 

Founc 

ation  Characteristics 

Winkler 

type 

Continuous 
type 

Mixed  mode 
(natural  soil) 

2 

1 

SxlO*^ 

1 

k 

s 

inland 

1.6 
1.2 

2.25 
22.5 

3x10^ 
6x10^ 

50 
100 

Ib/in^ 

lb 

soil  and 
sea  floor 

1. 

225. 

9x10^ 

500 

8.6 

18.3 

were 

.9 

1000 

500 

1000 

chosen 

.85 

1500 

as 

.8 

2000 
2500 

3000 

6000 

examples 

3000 
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Appendix  B  contained  the  tabulated  results  of  Figures  9  and  10  as 
well  as  data  on  the  two  foundations,    inland  soil  and  sea  sediment.      The 
results  show  that  the  difference  between  the  Winkler  hypothesis  and 
other  hypotheses  is  negligible  for  simply  supported  beams  with  uniform 
properties  and  loads.     This  may  not  be  true  for  beams  with  other 
boundary  conditions. 

B.      A  FREE-FREE  NON-UNIFORM  BEAM  PROBLEM 

To  show  how  the  Finite  Element  Program  can  be  used  for  a  general 
case  we  consider  the  following  problem: 


lb/in 


> 


400 
lb/in- 


k  =  kj   +  k^  V 


50.    in 
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A  symmetric  problem  was  considered  only  to  minimize  the 
computer  effort.      The  result  is  given  by  Table  5  and  Figure  11. 
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Table  5.     Deflection  of  a 
Free-Free  Non-Uniform  Beam 


£ 


I 

Location  #1  2  3  4  5 

I 

.  _  _  -  - 


Location  #  Deflection,    inches 

1  1. 19150 

2  2.09179 

3  2.87908 

4  3.40052 

5  3.58134 
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Deflection, 
inches 


0. 


1. 


2. 


4. 


Location  # 


Figure  11.      Deflection  of  a  free-free 
non-uniform  beam 
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VII.     CONCLUSIONS  AND  RECOMMENDATIONS 

The  computer  program  presented  in  this  thesis  provides  an  accurate 
and  reliable  means  for  solving  a  variety  of  one-dimensional  problems  of 
beams  on  nonlinear  foundations.     The  use  of  this  program  and  efforts  to 
increase  its  versatility  are  encouraged. 

A.      REMARKS  ON  THE  PROGRAM 

The  following  observations  have  been  made  during  the  development 
of  the  program  and  the  investigation  of  various  problems 

1.       Curve  Fitting 

It  was  assumed  that  the  reaction  of  the  foundation  is  proportional 
to  V*^      ,    where  V  is  the  deflection  and   p     any  positive  real  number. 
Since  the  Finite  Element  Method  via  the  Galerkin  procedure  is  not  easily 

applied  directly  for  the  case  when    p     is  not  an  interger,    V^        has  been 

2 
replaced  by  a  least  square  fit  polynomial  of  the  form  b     +  boV  +  b   V    . 

The  reaction  of  the  foundation  must  be  zero  when  there  is  no  deflection, 

therefore  the  greatest  error  may  come  from  the  non-zero  residual 

constant  b     (Figure  12).     Suppose  the  least  square  fit  curve  r^  =  b,   + 

b_V  +  b   V       intercepts  the  curve  r     =  V  at  the  point  V  =  V.     If  the 

deflection  of  the  foundation  falls  in  the  vicinity  of  V  ,    the  result  is  highly 

accurate.     The  result  is  less  accurate,    if  the  deflection  is  much  larger 

or  smaller  than  V  .     One  must  estimate  how  large  or  small  the  deflection 

will  be,    in  order  to  select  a  "best  fit"  curve. 
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r     =  kV 
n 


rf  =  bi+b2V+b3V 


V 


Figure  12.     Accuracy  of  curve  fitting 

2.  Initial  Estimate  of  Deflection 

Iterative  procedures  are  among  the  more  popular  schemes  for 
solving  systems  of  nonlinear  algebraic  equations.     Such  methods  require 
an  initial  estimate  of  the  solution.     In  this  work  the  initial  estimate  was 
taken  in  the  form  of  a  second  order  polynomial  in  the  subroutine 
VGUESS.     Not  all  initial  estimates  will  yield  a  solution;  therefore  it  is 
desirable  that  if  one  initial  estimate  fails,    it  be  replaced  by  another 
one.      This  is  accomplished  in  subroutine  VGUESS  which  changes  the 
scale  of  the  initial  estimate  function. 

3.  Computational  Considerations 

In  any  computer  effort  there  are  two  major  concerns,    the 
space  required  for  a  program  and  the  time  necessary  for  solution. 
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For  the  linear  Finite  Element  problem  the  storage  area  required  is 
proportional  to  m  X  M,    where  m  is  the  number  of  coordinates  and  M 
is  the  bandwidth  of  the  system.     A  change  in  m  or  M  yields  a  correspond- 
ing linear  change  in  the  size  of  the  computer  storage.     In  addition  the 

time  for  solution  (for  a  Gauss  elimination  type  scheme)  is  of  the  order 

2 

of  m  X  M     for  a  symmetric  system.      These  considerations  show  the 

disadvantages  which  result  when  m  or  M  increase,    and  are  well  known. 

In  the  case  of  nonlinear  problems  the  effects  of  increasing  m 
and/or  M  on  computational  effort  are  even  more  pronounced.      For 

example  the  space  requirement  associated  with  the  nonlinear  terms  is 

2 

m     X  M.     Hence  an  increase  in  m  from  m     to  m     increases  the  storage 

in  the  ratio    { zLS    .      This  is times  greater  than  in  the  linear  case. 

mj^  rn. 

Moreover,    the  computational  effort  is  greatly  affected  since  a  large 
number  of  terms  greatly  increases  the  number  of  iterations.      For 
example  in  the  solution  of  an  8-element  model  (i.  e.  ,   m  =  18,    M  =  6) 
non  linear  problem,    the  storage  for  the  array  is  about  8000  bits  (single 
precision)  and  the  computer  time  was  about  60  seconds.     When  the  same 
problem  was  modeled  by  a  4-element  model  (i.  e.  ,    m  =   10,    M  =  6),    the 
storage  was  about  2400  bits  (single  precision)  and  the  computer  time 
was  about  6  seconds;  a  reduction  in  storage  in  the  ratio  of  3  to  1  and  a 
reduction  in  computational  tinne  in  the  ratio  of  10  to  1. 
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B.      REMARKS  ON  THE  WINKLER  HYPOTHESIS 

Let  r     =  kV^       be  the  total  settlement  curve,    i.  e.    the  deflection 
n 

versus  reaction  curve  of  a  real  foundation,    which  can  be  obtained  by  an 

experiment  [Ref.    30].      Let  the  straight  line      n  =  kV,    where  k  is  the 

foundation  modulus,    represent  the  behavior  of  the  Winkler  foundation 

(Figure  13).     Let  V    correspond  to  the  point  at  which  the  two  curves  r 

and  r     intersect.      Consider  the  area  under  the  curve  r     =  kV*^ 
n  n 

between  0  and  Vq.      This  area  is  equal  to: 


rV 


0 

rdV 


0 
Consider  the  functional 

•L 


/V 


0  -n  b-  P+1 

kvP    dv    =  _JL_  v^^  (3o; 

p+1      ^ 

0 


F=      I      |(EIV)        +  (S   V)       +JS V^"*"^     -  qV)ldx  (31) 

Jo    ^  2  P+1  f 

where  L  is  the  length  of  the  beam.      The  general  field  equation  (29)  is 

the  Euler's  equation  of  the  above  functional.      Therefore,    to  solve  the 

equation  (29)  is  the  same  as  to  minimize  the  functional  (31).      This 

functional  represents  the  total  potential  energy  of  the  system,    including 

the  beam,    the  foundation  and  the  applied  load. 

The  third  term  under  the  integral  is  the  energy  due  to  the  reaction 

of  the  foundation.     According  to  Eqn.    30,    the  energy  (31)  is  a  function 

of  the  area  under  either  the  r    or  r     curve,    depending  on  which  founda- 

1  n 

tion  is  being  used.     Neglecting  the  effect  due  to  the  shear  force  in  the 
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r 
lbs /in 


Winkler  foundation 


nonlinear  foundation 


V,    inches 


Figure  13.     Comparison  of  a  Winkler  foundation 
to  other  types  of  foundations 


foundation  for  a  while,    for  a  specified  beam  and  a  given  applied  load, 
one  can  see  that  if  the  deflection  is  small,    say  less  than  V,    the  Winkler 
foundation  deflects  more  than  the  nonlinear  foundation  does  for  equal 
supporting  effort.     Now,    if  the  effect  of  the  shear  force  in  the  foundation. 
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i.  e.  ,  —  V     ,    is  taken  into  account,    one  can  expect  that  the  deflection 

'2 
should  be  less.      The  shear  force  of  the  foundation  depends  on  (sV) 

"2 
which  is  usually  much  smaller  than  (EIV)        and  may  be  neglected.      This 

discussion  agrees  with  the  results  in  two  examples  of  two  different 

types  of  soils  given  in  Chapter  VI,    part  A. 


C.      RECOMMENDED  FUTURE  STUDIES 

1.  Following  a  similar  procedure,    the  problem  of  plates  or  shells 
on  two-  or  three-  dimensional  nonlinear  foundations  may  be  investigated. 

2.  This  thesis  dealt  with  nonlinear  foundations  which  are  either 
continuous  or  discontinuous.     Real  foundations  may  be  some  combination 
of  these  foundations;  hence  there  is  a  need  to  bridge  the  gap  and  consider 
foundations  where  the  deformation  is  partly  localized  and  partly 
continuous. 

3.  In  the  present  investigation  the  nonlinear  term  kV  was 

2 
approximated  by  the  second  order  polynomial,    b,    +  b2V  +  b    V    .      This 

restriction  results  from  the  finite  element  approximation  of  V(x)  by 

m 

N        V.G.(x).      Future  studies  might  consider  ways  to  treat  kV^         or  in 

1=1 

fact  the  more  general  case  of  an  arbitrary  nonlinear  function  of  V, 

directly. 
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APPENDIX  A 

CONTINUOUS  FOUNDATION  ANALYSIS 

A.      GENERAL  FORMULATION 

This  formulation  is  limited  to  plane  stress  theory  only.     Let  u(x,  y) 

and  V  (x,  y)  be  the  displacements  at  a  certain  point  of  the  foundation,    in 

the  X-  and  y-  directions.     In  general,    an  approximate  solution  can  be 

obtained  by  expanding  u(x,  y)  and  v(x,  y)  in  a  finite  series: 

m 

u(x,y)  =22    Ui(x)     <p.(y)    , 

(32) 


i=l 
n 


v(x,y)  =Z^    V  (x)      ^.(y)      , 
j  =  l  J 


where  the  dimensionless  functions        <p.(y)  and        ^(y)  are  known  and 
the  functions  U.(x)  and  V-(x),    which  have  the  same  dimensions  as 
u(x,  y)  and  v(x,  y),    are  unknown. 

From  the  theory  of  elasticity,    in  the  two-dimensional  plane  stress 
case,    the  stresses  and  strains  are  related  as  follows  [Ref.    3]: 


'XX 


'yy 


i..j 


xy 


1-^ 


yx 


<     ^xx      ^  n^YY-    ^ 


(     u^ 


f      ^xx 


E. 


'      ^YY    ' 


2(l+u^) 


xy 


(33) 
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XX 


•yy 


'xy 


dn 

ay 

"^yx     = 


> 


(34) 


b^     j_      av 


ay 


ax 


where  E ,.  and    ^r  are  the  modulus  of  elasticity  and  the  Poisson's   ratio  of 

the  foundation,     e        and    ^       are  the  strains  in  the  x-  and  y-  direction 
XX  yy 

and   ip       is  the  shearing  strain. 

In  terms  of  the  series  (32),    the  relations  (33)  may  be  rewritten: 


m 


n 


'XX      = 


1-  P. 


irl 


Ty 


1-v^      trf  j-'       ^ 


> 


(35) 


yy\ 


'xy    = 


'yx 


■fZ^i  ^'i  ^  I/j  ^  3 


2(1+  l^f)  L--1  j--l 

where  the  "prime"  denotes  the  derivative  of  the  function  with  respect 
to  its  own  argument. 

From  the  foundation  model  shown  in  Figure  I,    consider  now  an 
elementary  strip  of  length       Ax  (Figure  14).      The  functions  U.(x)  and 
V.(x)  can  be  obtained  from  the  equilibrium  conditions.     According  to 
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V^y^T^^"" 


Zt"+  ^^AX 


XX"  ax 


Figure  14.      Force  system  on  a  strip 
of  the  foundation 
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Lagrange's  principle  of  virtual  displacements,      the  total  work  of  all 
internal  and  external  forces  acting  on  this  strip  due  to  any  virtual  dis. 
placement  is  zero.     Let: 

u  (x,y)  =  ^    ^    <P  i(y)      .  (36) 

V(x,y)  =    >^    \r       ^    (y)      ,  (37) 


5^  *i 


be  the  virtual  displacements,    where  U.  and  V-  are  arbitrary  constants 
at  that  location    x    of  the  strip. 

Since  there  is  no  force  applied  on  the  surface  of  the  foundation, 
and  the  bottom  of  the  foundation  rests  on  a  rigid  base,    there  is  no  work 
done  on  these  two  faces.     It  is  assumed  that  all  properties  of  the  founda. 
tion  remain  constant  in  the  z-direction,    hence  there  is  no  pressure 
gradient  in  the  z-direction  and  therefore  the  total  work  done  on  two 
faces  perpendicular  to  the  z-direction  is  zero.      The  system  of  forces 
which  must  be  taken  into  account  is  shown  in  Figure  14. 

The  external  forces  are  the  result  of  normal  stresses        q-     , 

OUvV  11- 

Q.  +    ■       Ax,    the  shearing  stresses      'T'xv'         ^xv 

XX  9  x  ,  y  / 

^'^xy        Ax  and  the  distributed  forces,    X  (x,  y)  in  the  x-direction  and 

Y(x,  y)  in  the  Y-direction. 

The  internal  forces  are  the  result  of  normal  stresses       n-       s^nd 

o-yy 

shearing  stresses       fj'     .      The  work  done  by  these  internal  forces  is 


A  virtual  displacement  is  any  displacement  which  is  consistent 
with  the  forces  and  constraints  [Ref.    2]. 
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the  strain  energy  of  the  strip  [Refs.    4,  9]. 

The  total  work  done  on  the  strip  by  any  m  +  n  virtual  displacements 

is  then  given  by  the  following  expression; 

H  H 

r       (    _^^  Ax)  (wdy)      u.     -  /*(   r^y  Yxy    )C^Kwdy) 


0 


0 


X    (x,y)(Axwdy)    u-     +  (^^^  A  x  )  (wdy)    Vj  (38) 


H  H 


r 

^  ^yy    ^yy^  ^  ^^   ^^^  ^  "^        J    Y  (X'  y)    (     ^x   wdy  )      v .     =    o 


0  0 

(i  =  1,2,  3,  .  .  .  ,  m      ;      j  =  1,  2,  3,  .  .  .  ,  n) 

where 


Y  +      L_      =        U.       <P      (y)  (39) 

»xv  ":\  -^  1 


xy 


■t)  V, 


f 


e  ...  L         =  ^  ,        ^      iv)  (40) 


yy  -.  J 

7>y 


f 


2)'^> 


^  t  XX  *>7^ 

( Ax)  (wdy)       u^  =    work  of  external  forces  (_l_i2lAx)(wdy); 

"^  X  "^x 

H 

I      ^        ^^      'xv^  (  ^^    wdy)  =  strain  energy  due  to  shearing  stress        Txy' 
Jo                      ^  -  ^ 


f 


■         -  -  -       -  -    -  ■B'^xy         Ax  )    (wdy)  ; 


( Ax)  (wdy)      ^-     =    work  of  external  forces     ( ^     xy 

'o         't^   X  "b  X 


(  £         )    (  Ax    wdy)     =    strain  energy  due  to  normal  stress 


f 

-'o      "^yy   ^yy  "yy 

I     X  (x,yX     Ax  wdy)  u.   +       /  i       r  i.     i      /• 

J  ^  I      Y(x,  y)(    Ax  wdy)    Vj      =  work  of  body  forces . 


After  dividing  through  oy  wAx  and  substituting  U. ,    v.  ,       y     ,     T^y  by 
(36),  (37),    (39)  and  (40),    the  expression  (38)  maybe  rewritten  as  follows: 


m   r  /-H 


S 


T>   cr. 


XX 


■^i    ^i(y)  dy 


+ 


)    Txy  U.       Cf>  iiy)  dy  +     J 


(y)  dy  +     JqX  (x,y)  U'i(y)dy 


EI  ^5 


^V        v.   iL..(y)dy 


I. 


C       Y     (x,y)  V.      vj/ j  (y)  dyl-  0 


H 


CTyy^j    f  j(y)   dy) 


(41) 


or 


A 


y^l        U.     [         1  •7)gxx    y.   (y)  dy  -     \      Txy(p(  (y)  dy  +       \       X  (x,  y)y>.  (y)dy] 

i=l 


+     Z-rf    (  "^j  t    Jq     '5-^   4^  j^^'^  ^y  -       -'o    ^y^  H'j   <y)^y  +    Jo   Y(x,y)f  .(y)dy] 


(42) 


Since  U.     and  V-  are  arbitrary,    the  above  expression  leads  to  the  two 
1  J  -' 

equations  in  the  x-  and  y-directions : 
\         '^CTxx     cp-(y)dy  -        \       Txy(^.(y)dy+  C       X  (x,  y)  c^.  (y)dy  =  0  (43) 

r  ^Tx 


^U.j(y)dy-  \        CTyy      v|;.(y)dy+        J         Y(x,  y)  vp  j(y)dy  =  0  (44) 
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Substitution  of  (35)  into  (41)  and  (42)  leads  to  a  system  of  m  +  n 
second  order  differential  equations  in  U^(x)  and  V.(x)  as  follows: 
j^H  ^  rn  m  ^ 

r^  E       ^      I   ^  ^  ,  fH 

_\        !_       [^Uifi      +XI!vj  4^j  ]<^k    dy    +  X      <pj^dy  =  0y45) 

Jo    2(1  +  V.)         i=l  j  =  l  -^0 


*H  -c^  rn  n         I 


'0 


"()x     2( 


1  +  ^f)        i-l.  j=l 


•"  *"  ,         n        ,  r^ 

i     -^    [   V   Zu>      +ZZ^if'     ]    4^^    dy    +  Y^     dy    =    0 

i^i-Vf  i=i    '  '       j-1    ■"    -^  *^0 

(k=l,2,3,...m 
i  =  1,2.3, ....n) 
These  are  the  two  governing  equations  of  a  two-dimensional 
foundation  which  has  a  depth  H. 

B.      SPECIAL  CASE 

For  the  limited  scope  of  this  thesis,    let  the  horizontal  displace- 
ments u  be  negligible.      This  means  that  any  virtual  displacements  in 
the  x-direction  are  also  negligible.     Since  U.     are  arbitrary  constants, 
Eqn.    36  implies  that     <P:(y)  must  be  identically  zero.     Only  the  last 
equation  of  the  system  (45)  remains  and  it  reduces  to: 
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Z{l^a^,i    ""^^^^^'^-c^     J^^.H'.r.^vJ 


Y(x,y)      vl/  .  dy  =  0  (46) 

0  '    ^ 


i=l 

-  •  i" 

where  j  =  1,  2,  3,  ,  .  .  ,  n,    and  Er  and  V.  are  assumed  to  be  constant. 
The  above  expression  may  be  rewritten: 

22'   \     -^  ^^        ^f  \       yitJ    dy)      v|    ]      _  [    Ef  ^^'.    ^/.    dy]    V. 

~^   \     '2)x         2(1+Vf)         Jo  l^      -^0 


+  I  Y  (x,y)  vp  .    dy    =    0 


"j- 


is  the  reaction  of  the  foundation. 

H 


k;:      =  Ef 


2(1  +  Vf) 

H 


^0 


ij 


l-V^ 
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(47) 


or 

m 

-5  =  IZ      {(^ijvi)'     -    kijVij  (48) 

j  =  1,  2, 3,  .  .  .  ,  n 
where 


\  Y  (x,y)  c^  .  dy  (49) 

^0  ^ 


^ij  =  ^- )       Yi     Yj        dy  (50) 


dy  (51) 


The  coefficients  S.^,    lbs/in,    and  k--,    lbs/in   ,    are  the  two  characteristics 
of  a  foundation  [Ref.    5].      The  former  accounts  for  the  shearing  stress 
distribution  in  the  foundation  material,    and  the  latter  is  equivalent  to  the 
Wihkler  constant  of  proportionality.      They  are   however,    dependent  on 
how  the       jp.  functions  are  chosen  and  are  consequently  not  independent 
of  each  other. 

For  simplicity,  assume  that  the  foundation  has  an  infinite  depth, 
which  is  usually  true  for  any  soil  layer  compared  to  the  deflection  of 
most  structural  systems.     Let  the  one-dimensional  function  ^(y) 

satisfy  the  boundary  conditions: 

^(o)  =  1.  (52) 

4>{oo)  =  0.  (53) 

Then  the  displacement  at  any  point  of  the  foundation  is: 

v(x,y)  =  V  (x)     4>{y)  (54) 

in  which  V(x)  is  the  deflection  at  the  surface  of  the  foundation,    where 
y  =  0.      Furthermore,    assume  that  the  displacements  decrease  ex- 
ponentially with  the  depth  of  the  foundation,    then      ^^(y)  may  be 

expressed  as: 

-   Xy 
ipiy)  =  e      /  (55) 

where    X    is  a  known  constant  which  depends  on  the  properties  of  the 

foundation.     The  values  of    s    and    k*  are  then  given  by  (50)  and  (51): 

S    =  ^'  r     e-     "^dy=         ^i  lbs/i„    (56) 

2(1  +  Vf)  4X(1  +  Vf) 

*^0 
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k*=    _^f_         r°°      (e-^y)'^dy=:  ^^£  lbs/.    3      (57) 


^     Jo 


2(1.  V/  ) 


in 


Equation  (48)  may  be  rewritten 

r  =  (sV)'   -  kfV  (58) 

where  s  and  k*are  given  by  (56)  and  (57). 

C.      NUMERICAL  APPLICATIONS 

In  accordance  with  soil  mechanics,    taking  the  length  of  a  loaded 
area  as  unity,    Reference  7  gives  an  empirical  relation  between  bearing 
load  and  settlement  of  a  foundation  as  follows: 

r  =  1     L-    V  (59) 

^      1-Vf 

where    ^     is  called  the  influence  coefficient,    which  depends  on  the 

shape  of  loaded  area  and  the  properties  of  the  foundation.      The  values 

of       pt       are  given  in  Table  6,  [Ref.    6],     Comparison  of  the  Winkler 

hypothesis 

r  =  kV  (60) 

to  equation  (59)  yields: 

k=    _±         ^f  (61) 

M       1-V^f  -        - 

Eqns.    57  and  61  show  the  similarity  of  the  theoretical  analysis  and 

the  empirical  result.      The  value       X     can  be  computed  by  equating 

(57)  to  (61): 

X  =  2M  (62) 


71 


Table  6 
Influence  Coefficients 

(After  Z.    Wilun  and  K.    Starzewski  [Ref.    7]) 


Shape  of 

Flexible  foundation 

Rigid  foundation 

Settlement 

Settlement 

Mean 

foundation 

of  centre 

of  corner 

Settlement 

Settlement 

A^O 

^c 

M 

Mr 

0.  64  point  on 

Circular 

1.00 

perimeter 

0.85 

0.79 

Square 

1.  12 

0.95 

0.88 

Rectangle 

Mc^Mo 

2  X  1 

1.53 

1.30 

1.22 

5  X  1 

2.  10 

1.83 

1.72 

10  X  1 

2.53 

2.25 

2.  12 

100  X  1 

4.00 

3.69 

-- 
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For  the  purpose  of  this  thesis,    consider  a  loaded  area  of  shape 
wx      A  X,    where         Ax  is  taken  as  a  unit  length,    one  inch,    and  w  is 
equal  to  the  width  of  the  foundation  (Figure  14).     Assume  that  w 
ranges  between  5  inches  and  10  inches.      From  Table  6,    1.  83  <  ^  <  2.  25. 
Let    y.     =  2,    then,    from  (62)         \    =  4.     in  order  to  determine     s    and    k, 
the  modulus  of  elasticity  and  the  Poisson's  ratio  of  the  foundation  must 
be  known. 

The  modulus  of  elasticity  of  rocks  ranges  frora  100    to  10,  000 

2  2  2 

mn/m   (7mn/m     =     Iklbf /in   )  and  the  Poisson's  ratio  from  0.  1  to  0.  3 

[Refs.    10,    7].     Deep  sea  sediment  cores  have  modulii  of  elasticity 

ranging  from  0  to  12.4  psi  [Ref.    11].     The  modulii  of  elasticity  of 

beach  soils  vary  from  30  to  570  psi,    while  those  of  inland  soils  vary 

from  3,000  psi  to  11,500  psi  [Ref.    7].      Most  soils  have  a  Poisson's 

ratio  near  to  0.  5  [Ref.    12].      The  values  of    s    and    k    computed  from 

(56)  and  (57)  for  various  types  of  soils  are  given  in  Table  7. 

Table  7 
s    and    k    Values  Depending  on  Various  Types 


of  S, 

oils  (Loaded 

A 

rea 

of  Shape 

Ax:  = 

1". 

5"    < 

w 

<      10") 

Sea  Floor 

Beach  Soi 

Is 

General 

Sediment 

Inland  Soils 

Ef,    psi 

12.4 

30-570 

3000-11, 500 

^^f 

0.  5 

0.  5 

.45 

s,    lbs 

18.6 

50 

5,000-21, 000 

k,    psi 

8.3 

20-380 

1800-7200 
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APPENDIX  B 
DATA  AND  RESULTS  FOR  PROBLEMS  OF  CHAPTER  VI,    PART  A 

The  following  example  refers  to  Figure  23-10  of  Reference  30. 
Assume  that  the  piles  used  for  this  soil  have  24-inch  diameters.      The 
reaction  for  the  foundation  is 


r  =      (Test  load)  w  lbs /in 

ttD    /4 

With  D  =  24  inches  and  w  =  6.  5  inches  we  have  the  following  table. 


Table  8 
Deflection  Versus  Reaction  of  a  Natural  Soil 

From  Experimental  Data 

[Ref.    30]  lbs /in 

Deflection  Test  load 

(inches)  (tons) 

0.10  20.  574.6 

0.25  40.  1149.2 

0.40  50.  1436.5 

0.60  60.  1724.5 

0.90  70.  2011.8 

1.30  73.  2097.6 

A  curve  of  this  data  is  shown  in  Figure  15. 
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r(kip/in) 

r(lbs/in) 

inland 

sea  floor 

soil 

model 
9. 

3. 

/ 

• 

/ 

^^^■^ 

6. 

2. 

^ 

--Q 

_ — — — 6 

1.. 

> 

3. 

c 

0. 

1/ 

( 

). 

.5 

f       « ■  — 

1. 

V(in) 

(j_y  Experimental  data  curve 
M  Least  square  fit 

\J  Assumed  foundation  behavior  curve 

Z\  Approximate  Winkler  hypothesis,    k  —   3000.  lbs /in 

Figure   15.      Foundation  reaction  versus  deflection 
(After  Spangler,    M.  G.  ,    [Ref.    30]) 
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A.  A  TYPICAL  INLAND  SOIL 

Figure  15  shows  reactions  versus  displacements  of  a  typical  inland 

soil  (read  scale  on  left  hand  side). 

2 
Assume  that  the  Winkler  foundation  modulus  is  about  3000  lb/in    . 

The  deflection  results,    obtained  by  the  Finite  Element  program,    of 

a  simply  supported  uniform  beam  with  uniform  load  resting  on  this  type 

of  soil  are  given  by  Table  9. 

Column  1:    Solution  of  (EIV")"  +  kV  =  q 

II 
Column  2:    Solution  of  (EIV")     +        r{V)  =  q  where    r    (V)  is  replaced 

by  a  least  square  fit  curve  of  the  reaction  versus  deflection  curve  given 

in  Figure  15. 

ti  ,   I 

Column  3:     Solution  of  (EIV")     -  (sV  )     +  r{V)  =  q  where    s    and    k 

corresponding  to  a  particular  foundation  is  given  by  Table  7. 

B.  AN  ASSUMED  SEA  FLOOR  MODEL 

Since  there  is  insufficient  experimental  data  of  sea  sediment,    it 
will  be  assumed  that  the  reaction  versus  deflection  curve  of  sea 
sediment  has  the  same  shape  as  that  of  inland  soil.     According  to 

Table  4  and  Reference   11,    assume  that  the  Winkler  foundation  modulus 

2 
of  sea  floor  is  about  8.  3  lb/in    .     It  will  be  assumed  that  the  inland  soil 

reaction  versus  deflection  curve  is  scaled  down  by  8.  3/3000  to  obtain 

the  sea  floor  reaction  versus  deflection  curve.      (Figure   15.      read  right 

hand  side  scale). 
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The  deflection  results  of  a  simply  supported  uniform  beam  with 
uniform  load  resting  on  this  type  of  sea  foundation  model  is  given  by 
Table  10. 
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Table  9.     Compared  Foundation  Deflection 
of  Various  Hypotheses  for  a  Typical  Inland  Soil 
Represented  in  Figure  15 

Q  2 

a)       L  =  100  in,    EI  =  3  x  10     Ibs-in   ,    q  =  22.  5  lbs/in 

k  =  3000  lbs/in^,    s  =  6000  lbs/in 


Location 

(2)  Modified 

(3) 

1  Mixed  mode 

x/L 

(1) 

1  Winkler 

Winkler 

foundation 

0.0 

0.0 

0.0 

0.0 

0.  125 

0.003690 

0.350190D. 

.02 

0.  349588D-02 

0.250 

0.006390 

0.603699D- 

.02 

0.  602663D-02 

0.325 

0.007898 

0.744817D. 

.02 

0.743554D-02 

0.500 

0.008370 

0.  788845D. 

.02 

0.787515D-02 

7  2 

b)       L  =  100  in,    EI  =  3  xlO    Ibs-in   ,    q  =  22.  5  lbs/in 

2 
k  =  3000  lbs/in   ,    s  =  6000  lbs/in 


Location 

(2) 

1  Modified 

(3) 

1  Mixed  mode 

x/L 

(1) 

1  Winkler 

Winkler 

foundation 

0.0 

0.0 

0.0 

0.0 

0.  125 

0.005513 

0.524425D. 

.02 

0.  522121D-02 

0.250 

0.007740 

0.728582D- 

.02 

0.  726411D-02 

0.375 

0.007988 

0.749146D. 

.02 

0. 748313D-02 

0.500 

0.008010 

0.739400D. 

.02 

0. 739234D-02 

7  2 

c)        L  =  100  in,    EI  =  3  X  10     Ibs-in   ,    q  =  225  lbs/in 

k  =  3000  lbs/in    ,    s  =  6000  lbs/in 


Location 

(2)  Modified 

- 

(3)  Mixed  mode 

x/L 

(1)  Winkler 

Winkler 

foundation 

0.0 

0.0 

0.0 

0.0 

0.  125 

0.055125 

0.536074D- 

.01 

0.533658D-01 

0.250 

0.077400 

0.749189D- 

.01 

0.746865D-01 

0.375 

0.079875 

0.773826D. 

.01 

0. 772876D-01 

0.500 

0.080100 

0.764925D- 

-01 

0. 764673D-01 
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Table  10.      Compared  Foundation  Deflection 

of  Various  Hypotheses  for  a  Sea  Floor  Soil 

Represented  by  Figure   15 

7  2 

L    =  100  in,    EI  =  3  X  10     Ibs-in 

a)  q  =  1  lbs /in 


o< 

;ation 

X 

:/L 

0. 

0 

0. 

125 

0. 

250 

0. 

375 

0. 

500 

Winkler 

0.0 

0. 131617E-01 

0.241054E-01 

0.312727E-01 

0.337603E-01 


Modified 
Winkler 

0.0 

0.  130091D-01 

0.238239D-01 

0. 309054D-01 

0. 333630D-01 


Mixed  mode 
foundation 

0.0 

0.  130029D-01 
0.238124D-01 
0. 308904D-01 
0.333467D-01 


b)  q  =  2.25  lbs/in 


Location 

X 

7L 

0. 

0 

0. 

125 

0. 

250 

0. 

375 

0. 

500 

Winkler 


0.0 

0.296138E-01 

0.542371E-01 

0.703635E-01 

0.759606E-01 


Modified 
Winkler 


0.0 

0.293785D-01 
0. 538043D-01 
0.698002D-01 
0.753519D-01 


Mixed  mode 
foundation 


0.0 

0.293644D-01 
0.  537781D-01 
0.697660D-01 
0.753150D-01 


c)  q  =  4.  5  lbs/in 


Location 

Modified 

Mixed  mode 

x/L 

Winkler 

Winkler 

foundation 

0.0 

0.0 

0.0 

0.0 

0.  125 

0.592276E-01 

0.  591543D-01 

0.591253D-01 

0.250 

0.  108474 

0.  108346 

0. 108292 

0.375 

0. 140727 

0. 140568 

0. 140498 
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d)  q  =  6.75  lbs/in 


location 

Modified 

Mixed  mode 

x/L 

Winkler 

Winkler 

foundation 

D.O 

0.0 

0.0 

0.0 

0.  125 

0.8884616E. 

-01 

0.893441D. 

.01 

0.  892998D-01 

0.250 

0. 162712 

0.  163656 

0.  163574 

0.375 

0.211090 

0.212344 

0. 212237 

0.500 

0.227882 

0.229248 

0.229132 

e)  q  =  9  lbs /in 


Location 

Modified 

x/L 

Winkler 

Winkler 

0.0 

0.0 

0.0 

0.  125 

0.  118455 

0. 119966 

0.250 

0.216949 

0.219677 

0.375 

0.281454 

0.285174 

0.500 

0.303843 

0.307884 

f)  q  =  15.75  lbs/in 


Location 

Modified 

x/L 

Winkler 

Winkler 

0.0 

0.0 

0.0 

0.  125 

0.207297 

0.214630 

0.250 

0. 379660 

0.393300 

0.  375 

0.492545 

0. 510477 

0.  500 

0.531724 

0.551182 

g)  q  =  22.  5  lbs/in 


Location 

Modified 

x/L 

Winkler 

Winkler 

0.0 

0.0 

0.0 

0.  125 

0.296138 

0.313991 

0.250 

0. 542371 

0. 375550 

0.375 

0.703635 

0.747221 

0.500 

0.759606 

0.806886 

Mixed  mode 
foundation 

0.0 

0. 119906 
0.219657 
0.285028 
0. 307726 


Mixed  mode 
foundation 

0.0 

0.214517 
0. 393091 
0. 510204  . 
0.550887 


Mixed  mode 
foundation 

0.0 

0.313816 
0. 575277 
0. 746800 
0. 806430 
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APPENDIX  C 

COMPUTER  PROGRAM  NOMENCLATURE 

A  complete  listing  and  description  of  all  variables  used  in  the 
program  is  not  practical.      The  items  listed  in  this  appendix  are  common 
to  several  areas  of  the  program  and  will  assist  the  reader  in  a  study  of 
the  program.      For  convenience  items  are  listed  by  variable  type.      The 
definition  or  function  and  dimension,    if  applicable,    of  each  item  is 
given. 

A.      INTEGER  CONSTANTS 

KA  maximum  number  of  data  of  foundation  modulus 

KEI  maximum  number  of  data  of  beam  flexural  rigidity 

KQ  maximum  number  of  loading  conditions 

KS  maximum  number  of  various  type  of  soils 

KATEST  specified  foundation  modulus  for  test  program  only 

KETEST  specified  beam  flexural  rigidity  for  test  program  only 

KQTEST  specified  loading  condition  for  test  program  only 

ITMAX  maximum  number  of  iterations 

NBC  number  of  boundary  conditions  imposed 

NCODE  code  number,    if  NCODE  is  not  equal  to  1,    the  shear  force 

in  the  foundation  material  is  not  taken  into  account 
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NDOF  number  of  degrees  of  freedom 

NELEM  number  of  elements 

NELMAX         maximum  number  of  elements 


NNOD 
NSIG 


number  of  nodal  points 
number  of  significant  figures 


B.      REAL  CONSTANTS 


foundation  modulus  k 

beam  flexural  rigidity 

applied  load 

shear  force  in  a  soil  foundation 

total  length  of  the  beam 

beam  element  length 


ALPHA 

EI 

Q 

SSOIL 

TL 

X 

C.  VECTORS 

ASAV(IO)  a  set  of  foundation  modulii 

B(10)  coefficients  of  a  least  square  fit  polynomial 

EF(IO)  element  force  vector 

EISAV(IO)  a  set  of  beam  flexural  rigidities 

IDBC(20)  identification  number  of  boundary  conditions 

PO{10)  a  set  of  power  p's  in  the  term  kV 

QSAV(IO)  a  set  of  applied  loads 

SK(IO)  a  set  of  foundation  modulii  of  natural  soils 

SS(IO)  a  set  of  shear  forces  of  natural  soils 
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SYSF(20) 
THETA(20) 
V(20) 
WA{600) 

D.      MATRICES 
DC(4,4) 
DM{4,  4) 
DN(4,  4,  4) 
EK(4,4) 

ICORR(20,4) 
SYSC(20,  20,20) 
SYSK(20,20) 


system  force  vector 
slopes  at  the  nodal  points 
deflections  at  the  nodal  points 
working  vector 


Cj^j^  -  matrix,    referred  to  Eqns.    23  and  24 
M-1^  -  matrix,    referred  to  Eqns.    23  and  24 
^iik  "  n^^atrix  referred  to  Eqns.    23  and  24 
element  stiffness  matrix  K.,    referred  to  Eqns.    23 
and  24 

correspondence  table  of  local  and  global  points 
system  N^jj^  matrix,    referred  to  Eqn.    24 
system  Cj^.,    M^j^  and  K-i  _  combined  matrix, 
referred  to  Eqn.    24. 
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//TAM    JOB    ( 1714,0733»NS34), •  •,TIME=1 

//    EXEC    FGRTCLb,  REGIOr4=200K 
//FJRT.SYSIN    DO    ^ 
C 

c 

C  THIS    lAAlN    PROGRAM    MAY    BE    USED    FOR    EITHER 

C  LINEAR    UR    NONLIiMcAR    PROBLEMS     ,     ACCCRDING    TO 

C  APFROPKIATE    INSTRUCTION    .     FOR    LINEAR 

C  PRCBLEMSiOMIT    ALL     'IMPLICIT    REAL-8 ( A-H, G-Z  )  ' 

C-  CARDS     in    THE    MAIN    AS    WELL    AS    IN    ANY    SUBROUTII 

C 

c 

IMPLICIT    REAL*8( A-H,0-Z) 

COMMON     EK(4t4)  ,DC(4,4)  ,DM('*.4J  ,DN(4t^»^)  » 

1  IC0RR{20,4)  ,SYSK( 20,20) »SYSC( 20,20 ,20)  t 

2  SYSF( 20) ,tF (4) , I0BC(20) ,B( lOJ , 

3  NNOD , HE  LE  M , NDGF , N BC , T L , AA 1 , AA2 , 

4  X,EI ,Q,ALPHAtSS0IL,A1 , A2 
DIMENSION  WA(600) ,THbTA(20) ,V( 20)  , 

1  PO(IO) ,SK(10)  tSS(lO)  , 

2  ASAV( 10) tEISAVI 10) ,QSAV(10) 
EXTERNAL    AUX 

1  FORMATS 'l* ) 
C 

C      THE  NCODE  NUMBER  IS  USED  TO  CONTROL  THE 
C      OUTPUT  RESULTS  DEPENDING  ON  THE  TYPES 
C      OF  PROBLEMS  .  NCCDE^l  FOR  LINEAR 
C      PROBLEMS  .  NC0DE=2  FOR  NONLINEAR  PROBLEMS. 
C 

NCODc-2 

CALL    STOREiNDDF,NBC,TL,NELEM, 

1  KATEST,KETESTtKQT£ST,KA,KE1, 

2  KQ,KS»EPS,NSIG, 

3  ASAV,QSAV,EISA\/,P0,SS,SK,  ICORR) 

c 

C      THE  FOLLOWING  PARAMETERS  MAY  BE  USED  IN 

C      THE  CASE  OF  UNIFORM  BEAM  ^ITH    UNIFORM 

C      LOAD  AND  MAY  VARY  ACCORDINGLY  TO  THE 

C      BEAM  FLEXUkAL  RIGIDITY  ,  LOADING  CONDITION 

C      AND  TYPES  CF  FOUNDATIONS, 

C 

IP=1 

JP  =  1 

KP  =  1 

LF=1 

PRINT  1 

CALL    INCHK(NELEM.NNOD,NBC,NDOF, 

1  X,TL,EI,0, IP, JP,KP,ISOIL, 

2  SSOIL, ALPHA, NCCOE,VSCALE, 

3  ASAV,EISAV,QSAV , SK t SS , I DBC , V , PO) 
C 

C      FOR  THE  CASE  OF  LINEAR  PROBLEMS , REPLAC E 
C      'CALL  STIFF2«  CARD  BY  'CALL  STIFF!'  CARD 
C 

CALL  STIFF2 

CALL  LOAD 

CALL  BOUND 
C 

C      FOR  THE  CASE  OF  LINEAR  PROBLEMS , REPLACE 
C      THE  FOLLOWING  CARDS  BY 
C      CALL  SOLVEiSYSK, SYSF, V,NNOO) 
C      CALL  R£SULT(N^GD,V,£I,  lER,  ITMAX, NCODE) 


C 


NG  =  0 

CALL  VGUESS(NNOD,NELEM,VSCALE, ITMAX.V) 

CALL  ZSYSTM( AJX,EPS,NSIG,NNOD, V, ITMAX,WA ,PAR,IER) 

IF(IEK.EQ.O)  GO  TO  5 

NG  =  NG  +  I 

IF(NG.LE. 10)  50  TO  4 

CALL    RbSUL  i(NNOi)  ,V,EI  ,  lER,  IT  MAX,  NCODE) 

STOP 

END 
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SUBROUTINE    STORE ( NDGF , N8C» TLtNELEM , 

1  KATEST,K£TEST,kQTEST,KA,KEI  , 

2  KQ,KS,EPS,.NISIG, 

3  ASAV,y5AV»EISAVtP0,SS,SK,ICCRR) 
C 

C      FOR  LINEAR  PRGBL EMStOM  I  T 

C      'IMPLICIT  KEAL^8( A-H,0-Z)«  CARD 

C 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION    ASAV(IO)  ,QSAV{10)  , £  I S AV(IO} t SS (1 0 ) , 
1  SK(  lOltPjdO)  ,IC0Rk(20,4) 

c 

C      THE  FOLLOWING  PARAMETERS  ARE  USED  FOR 
C      TESTING  THE  CONVERGENCE  OF  THE  PROGRAM. 
C 

NELMAX=8 

KETEST=l 

KQTEST=1 

KATEST=1 

EPS=1.0D~06 

NSIG=5 
C 

C      READ  IN  ALL  INFORMATION  NEEDED  TO  THE 
C      ENTIRE  FIELD  OF  THE  PROBLEM. 
C 

READ(5,lui  KA,KEI ,KQ,NELEM 

READ{5t3)  NO OF, NBC ,KPO,KS»TL 

3  F0PMAT(4I5»GI5.3) 
READ{5,4)     (PG( I ) , I=1,KP0J 
READ(5,4)     (ASAV( I ) ,1=1 ,KA) 
READi[5,4)     (QSAV(  I)  ,i=--l  ,KQ) 
READ(5,4)     (lISAv/(  I  ),I=1  ,KEI) 
READ{5,4)     (SK( I) ,I=1,KS) 
READ(5,4)     (SS( n .1=1 tKS) 

4  F0RMAT(4G1 5.3) 
C 

C  READ    IN    THE    CDRR ESPONDENCE       TABLE 

C  BETWEEN    LOCAL    AND    GLOBAL    COORDINATES 


C 


DO  11  I=1,NELMAX 
11  RE  AD (5, 10) iCORR{  I, 1 ),ICORR( 1,2), ICORRd  , 3) , ICORR( 1,41 
10  FORMAT (41 5 i 

RETURN 

END 
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SUBROUTINE  INCHK < NELEM ,NNOO , NBC » NDOF, 

1  X,TL,EI ,Q,IP, JP,KP, ISOIL, 

2  SSGlLtALPHAfNCOOE tVSCALF, 

3  ASAV,EISAV,QSAV,SK,SS, IDBC,V,PO) 

c 

C      FOR  LINEAR  PRGBLEMS,OM IT 

C      • IMPLICIT  REAL*8( A-H,0-Z) •  CARD 

C 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  ASAV(10),£ISAV(10)»QSAV(I0), 

1  ID6C{2O),V(2O),P0(10),SK(lO),SS(10) 
VSCALE=1.0-0o 

Q  =  Q5AV(  IP) 

EI=EISAVIJP) 

IF(NCUDE.NE.l)  GO  TO  7 
C 

C      THESE  TWO  FOLLOWING  STATEMENTS  IS  DEALING 
C      WITH  NATURAL  SOILS 
C 

SSCIL=SS< ISOILJ 

ALPHA=SK( ISOIL) 

GO  TO  9 
C 

C      SSCIL   MAY  HAVE  A  CERTAIN  VALUE  IF  THE 
C      SHEAR  FORCE  IN  THE  FOUNDATION  MATERIAL 
C      IS  TAKEN  INTO  ACCOUNT. 
C 

7  SSOiL=0.0 
ALPHA=ASAV(KP) 

9  NN0D=2-(NELEM+1) 
X=TL/(2*NELEM) 
PRINT  6 tNNUD,N£LEM,NBC,NOOF ,X 

8  F0RHAT(5X,'  NUM3ER  OF  NODAL  PO I NTS  =  ' , I  5 t / , / 
15X, 'NUMBER  OF  EL cMENTS= *  ,  I  5  ,/ ,  /  , 

2  5X,' NUMBER  OF  bOUNOARY  COND IT  I ONS= S I  5 , / , / , 
35Xt'NUMBER  OF  OEGREcS  OF  FREEDOM=«  ,  I  5 t / » / » 
45X,' ELEMENT  LENGTH^' t IX ,015 .6 , //// ) 

bcTA=(  ALPHA/ (4. -EI  ))*=<--.  ^5 
C 

C     THE  FOLLOWING  BOUNDARY  CONDITIONS  MAY  NOT 
C      BE  NECCESARY  IN  GENERAL  AND  CAN  BE 
C      READJUSTED  IN  THE  SUBROUTINE  BOUND. 
C 

IDBC(1J=1 

ICBC(2)=NNUD 

V(1)=0.0 

V(NNOD)=0.0 

RETURN 

END 
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SUBROUTINE  STIFF2 

IMPLICIT  REAL«8( A-H,0-Z) 
C 

C-  THIS  SUBROUTIME  IS  USED  TO  COMPUTE 

C  THE  ELEMENT  STIFFNESS  MATRICES  AND  THE 
C  SYSTEM  STIFFNESS  MATRIX  FOR  NONLINEAR 
C      PROBLEM  ,  ACCORDING  TO  THE  FOkMULA  {   ) 

C      ON  PAGE  

C 

COMMON  EK(4,4)  tOC(4,4),DM(4,4)  ,DN{4,'*,^)  , 

1  IC0RRU0,4)  ,SYSK(20,20),SYSC(2j,20t^0»» 

2  SYSF(20),EF(4), IDBC(2G J ,3(10)  , 

3  NNCOtNtLEMfNDOF  ,NBC , TL , AAl , AA2, 

4  X, EI, Q, ALPHA, SS0IL,AItA2 
DO  5  I=1,NNCD 

DO  5  J=l,NNOO 

SYSK( I , J)=0. JD+O 

DO  5  K=l,N.^lUL) 
5  SYSC(  It  J,KJ=O.OD*-0 

DO  6  N=l iNELEM 
C 

C      THE  FOLLOWING  FUNCTIONS   SFUNCT  ,  EFUNCT  ♦ 
C      62  ,  B3  ARE  DEALING  WITH  VARIABLE  SNEAK 
C      FORCE  IN  THE  FOUNDATION  MATE RI ALS , VARI ABLE 
C      BEAM  FLEXURAL  RIGIDITY  AND  NONLINEAR 
C      REACTION  OF  THE  FOUNDATION  .  THE  REACTION 
C      OF  THE  FOUNDATION  MAY  BE  WRITTEN  AS 
C      R(X)  =  B2-V  +  33-V--!==^2 
C 

SSOIL=SFUNCTiN,X) 

EI=EFUNCT(N,X) 

B(l)=O.ODfO 

B(2)=B2 (N,X) 

B(3)=B3(N,X) 

SS0IL=SSO[L/£I 

A1=aLPHA/EI 

AA1  =  A1-!'B(2} 
C 

C      THE  FOLLOWING  IS  THE  BEAM  ELEMENT 
C      STIFFNESS  MATRIX 
C 

EK(1,  l)=12.D+00*(  1. 0+00 /X=^* 3) 

EK(l,2)=6.D*-00*(  I.0+00/X--2) 

EK(  I  ,31  =-12.0  +  00=5=(  l.D<-0C/X':=*3) 

EK(l,4i=6.D<-00^(  l.D  +  00/X**2} 

EK(2f2) =H.0+00-( 1 .0+00/ X) 

EK(2,3) =-6.0+00* ( 1.0+00/X**2J 

EK(2,4J=2.Df00v(l.D+0  0/X) 

EK(3,3)=12.D  +  0  0-(l.D+00/X*=!'3) 

EK(3,4)=-6.0+0  0*(l.D+00/X*-2) 

EK(4,4i=4.0+00«(l.D+00/X) 

EK(2,1 J=tK(l,2) 

EK{:5,i)=EK(l,3) 

EK(3,2)=  eK(2,3» 

EK(4,  U=EK(1,4) 

EK(4,2)=£K(2,4) 

EK(4,3)=EKi3,4J 
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c 
c 
c 
c 
c 


c 
c 
c 
c 
c 


THE  FOLLOWING  IS  THE  ELEMEi^^T  FOUNDATION 
STIFFNESS  MATRICES  ,  ASSOCIATED  TO  THE 
LINEAR  TERMS 


DC(1 

OC(i 

DC  (I 

DCd 

DC(2,1 

DC (2,2 

DC(2,3 

DC(2,4 

DC{3,I 

DC(3,2 

DC(3,3 

DC(3,^ 

DC(4, 1 

DC(4,2 

0C(4,3 

DC(4,4 

DM(1, 

OM(I, 

DM(1, 

DMdt 

DM(2, 

DM  {2, 

DM(2t 

DM(3, 

DM(3,4 

DM(4,4 

DM(2,I 

DM(3,1 

0M(4,1 

DM{3,2 

0M(4»2 

DM (4, 3 


=1.20+00/ 
=0.10+00 
=-OC(i, 1) 
=0C(1,2 ) 
=DC(1,2) 
=0.133333 
=-DC(  it2) 
:-0. 03333 
=DC{1,3/ 
=  DC(2, 
=  OC(i 
=  DC(2 
=  DC(1 
=  0C(2 
=  0C(3 
=  DC(2 
=  .3714290 
=.0523810 
=.1285710 
=-.030952 
=.0095240 
=.0309520 
=  -.007143 
=  .3714290 
=-.052381 
=.0095240 
=DM(l ,2/ 
=DM(1,3] 
=  DMil,4) 
=  0,^(2,3) 
=0M(2t4i 
=DM(3,4) 


D+00«X 
3D+00*X 


) 
»l) 

,3) 
,4) 
,4J 
,4) 
,2) 


+oo*x 

+  00^''X*^2 
f  00-X 
0+0C*X**2 
+  00-X=f"!=3 
+  00*X=4^*2 
D+00*X**3 
+  00-X 
D+00*X**2 
<  00*X*^3 


THE  CORRESPONDING  COMPONENTS  OF  ALL 
STIFFNESS  MATRICES  ASSOCIATED  TO  THE  LINEAR 
TERMS  MAY  BE  AOOED  AS  FULLOWS 

DO  1  I=1,ND0F 

DO  1  J=1,ND0F 

DCd,  J)=SSOIL*DCi  I  ,J) 

DM(I , J) =AA1*0M(  i  t J) 

EK(I,  J)  =EKU  t  JJ+OCC  I,  J)+DM(  I,  J) 
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c 
c 
c 
c 


THE  FOLLOWING  ARE  THE  COMPONENTS  OF  THE 
MATRIX  ASSOCIATED  TO  THE  NONLINEAR  TERMS 


DN( 

1. 

1,  I) 

DN< 

2, 

1,1) 

ONI 

3, 

III) 

DN( 

4, 

1,1) 

DN( 

2, 

2,  I) 

DNI 

3i 

2,1) 

DN< 

4i 

2,1) 

DN( 

3i 

3,1) 

DN( 

4, 

3,  1) 

DN( 

4, 

4,  1  ) 

DN( 

2, 

2,  2) 

DN( 

3i 

2,2) 

DN( 

4, 

2,2) 

DN( 

3, 

3,2) 

ONI 

4, 

»3,2) 

DN( 

4i 

4,2) 

DN( 

3i 

3,  3  ) 

DN( 

4, 

r3,3) 

DN( 

4, 

4,3) 

DNI 

4, 

4,4) 

DNI 

1, 

2,  1) 

DNI 

li 

3,  1) 

0N( 

2i 

-3,1) 

DNI 

Ii 

p4,l) 

DN( 

2i 

4,  1  ) 

DNI 

3, 

i4,l) 

DNI 

1, 

1,2) 

DN( 

ii 

2,2) 

ONI 

1, 

-3,2) 

DNI 

ii 

-4,2) 

DNI 

2 

rl,2) 

DNI 

2, 

-3,2) 

ONI 

2, 

-4,2) 

ONI 

3 

-1,2) 

DNI 

3 

-4,2) 

ONI 

4, 

-1,2) 

DNI 

1 

rl,3) 

DNI 

i 

-2,3) 

ONI 

I 

,3,3) 

DNI 

1 

-4,3) 

ONI 

2, 

-  1,3) 

ON 

2 

»2,3) 

DNI 

2, 

-3,3) 

DNI 

[2, 

-4,3) 

ONI 

[3 

,1,3) 

DNI 

[3 

r2,3) 

ONI 

[3 

r4,3) 

DN 

14 

rl,3) 

ON 

[4 

,2,3) 

ON 

(1 

,1,4) 

DNI 

[  1 

-2,4) 

ON 

[I 

-3,4) 

ON 

(1 

,4,4) 

ON 

:2 

,  1,4) 

ON 

[2 

r2,4) 

ON 

(2 

,3,'*) 

ON 

(2 

,4,4) 

DN 

[3 

,  1,4) 

ON 

(3 

,2,4) 

ON 

t3 

,3,4) 

ON 

(3 

,4,4) 

ON 

(4 

,1,4) 

DN 

(4 

,2,4) 

DN 

(4 

,3t4) 

=.3071430+00*X 

=  .03849?D+00=:'X**2 

=  .0642360+00'^X 

=-.01/0610+00*X**2 

=  .  0065490  fOO'-'^X=!'*3 

=  .013a39D  +  00>i^X**2 

=  -.0035  71D  +  00=i'X*=.'=3 

=.064286DfOO*X 

=  -.Ol3889D  +  00=^X*«2 

=  ,003i750  +  OO^X=!'*3 

=  .001191D*-00*X^*4 

=  .0031750  +  00*X-i'*3 

=-.0007  94D+00*X**4 

=  .017a6lD+00*X=:^*2 

=  -.003571Df00'^X**3 

=  .0007940  +  jO-X'!'*4 

=  .j07143D  +  00=^X 

=  -.0384920  +  00=:=X=:^*2 

=.006349D+00*X**3 

=  -.001I91D*-00'!=X**4 

=  DM (2,1, 

-1  ) 

=  DM(3,1  , 

1) 

=DN(3,2, 

1) 

=uN(4,l , 

n 

=  Divi(4,2, 

1) 

=DN(4,3, 

1) 

=0N(2, 1 , 

1) 

=DN(2,2, 

rl) 

=  OiMi  3  ,  2  1 

1   ) 

=  UN  ( 4  ,  2  1 

1) 

=DN<2,2, 

1  ) 

=DNi3,2, 

.2) 

=0N(4,2, 

-2J 

=DN(3,2i 

,1) 

=  iJisl(4,3 

-2) 

=DN(4,2, 

-1) 

=  DN ( 3  , 1  , 

rl) 

=DN(3, 2i 

1) 

^ DN (3,3, 

-1) 

=DN(4,3, 

-1  ) 

=  ON ( 3 , 2  . 

-1) 

=DN(3,2, 

-2) 

=  0.M(3,i 

-2) 

=  DN  (  4 ,  3  , 

-2) 

=0N(3,3 

-1  ) 

=DN(3,3, 

-2) 

=DN(^,3, 

-3) 

=  Di\(  4,3 

rl) 

=DN(4,3 

-2) 

=DN(4,1 

-1) 

=  Dia4,2 

-1) 

=  DN(4,  3 

,1) 

=0N(4,4 

rl) 

=nN(4,2 

-1  ) 

= DNI 4, 2 

-2) 

=DN(4,3 

r2) 

=  UN ( 4  ,  4 

,2) 

^DN('t,3 

,1) 

=0N(4,3 

,2) 

=UN(4, h 

r3) 

=UN(4,4 

,3) 

=0N(4,4 

rl) 

=DN(4,4 

r2) 

=DN(4,4 

,3) 
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c 
c 
c 
c 
c 


THE    SYSTEM    STIFFNESS    MATRIX    IS    OBTAINED 
BY    COLLECTING    THE    ELEMENT    STIFFNESS 
MATRICES    OVER    ALL    ELEMENTS    OF    THE    SYSTEM 


J,K) 


AA2=A1*B(3) 

00    4    K=I,NOGF 

J=1,ND0F 

I=I,NDaF 

J,K)=Aa2*DN(I 

I=ltNL)OF 

ICGRRiNtI ) 

J=1,N00F 

ICORRlNt J) 
SYSK(NSlfNS2)=SYSK(NSI,NS2J 
DO    7    K=ltNDQF 
NS3=    ICORR(N,K) 

SYSC<N3l,NS2,NS3)=SYSC(NSl,NS2,NS3) 
CONTINUE 
RETURN 
END 


DO  4 
DO  4 
DN(I, 
DO  7 
NSl  = 
DO  7 
NS2= 


+    EK(I  ,Ji 


+    DN(I-,J,K) 


90 


SUBROUTINE   LOAD 
C 

C      FOR  THE  UINEAR  PROBLEMS  i  USING  SINGLE 
C      PRECISION, THE  'IMPLICIT  RE AL*8 ( A-H,0-Z ) • 
C      CARD  MUST  BE  OMITTED 
C 

IMPLICIT  aEAL*8( A-H,0-Z) 
C 

C      THIS  SUBROUTINE  IS  USED  TO  COMPUTE 
C      THE   ELEMENT   FORCE   VECTOR 
C      AND  THE  SfSTEM  FORCE  VECTOR 


C 


COMMON  EK(4,4) ,DC(4,4) ,0M(4,4) ,DN(4,4,^) , 

1  ICORR(20,4) ,SYSK(20,20) ,SySC(20,20,20), 

2  SYSP(20) ,EF(4),IDBC(20) »B(10) t 

3  NNODfNELEH, NOCF ,NBC , TL , AAi , AA2, 

4  X, EI, Q, ALPHA, SS0IL,AI,A2 
DO    1    I=l,NNOD 

.    SYSF(I)=0, 00+00 

DO    3    N=1,NELEM 

EF(1J=     .5Di-0  0*X 

EF<2)=     .O83333D+O0*X**2 

EF(3)=     .50+00-X 

EF(4)=    -.03"^333D+00vX**2 

Q  =  QFUNCT(N,X)-A.l  =5^8(1  ) 

EF(I)=Q^^EF{1) 

EF(2)=3-EF12) 

EF(3)=g-EF(3) 

EF(4)  =  Q^.^EF(4) 

DO    2    NT=1,N00F 
>    EF(NT)=EF{NT)/EFUNCT{N,X» 

I=ICORR(N, L) 

J=IC0RR(N,2I 

K=IC0RR{N,3) 

L=IC0RR(N.,4J 

SYSF( n  =SYSF(  I)f    EF( I) 

SYSF(  J)=SYSr(  J)     <■    EF(2) 

SYSF(K)=SYSF(K)     +    EF(3J 
I    SYSF(L)=SYSF(L)     <■    EF(4) 

RETURN 

END 
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SUBROUTINE    BOUND 
C 

C  FOR    LINEAR    PROBLEMS, THE    'IMPLICIT    REAL*8 

C  (A-H,0-ZJ'     CARD    MUST    BE    OMITTED 

C 

IMPLICIT    REAL*8( A-H,0-Z) 
C 

C  BEFORE     SOLVING    THE    ALGEBRAIC    EQUATION    SYSTEM 

C  THE    BOUNDARY    CONDITIONS    MUST    BE    APPLIED 

C 

COMMON    EK<4t^) ,DC(4,4) ,DM<4,4) ,DN(4,4,^) , 

1  ICORR(20,4}  ,SYS;<(  20,20)  ,SYSC(  20,20,  20 i  , 

2  SYSF(20),EF(4), I0BC(20) ,B(10) , 

3  NNOD.NdLEM, NDCF  ,NBC,TL, AAi, AA2, 

4  X, EI ,Q, ALPHA, SS0IL,A1,A2 
DO  3  N=1,NBC 

C 

C      FOR  FREE-FREE  ENDS  BEAM  IT  IS  DESIGNED 

C      TO  SKIP  THE  FIRST    BOUNDARY  CONDITION 

C      WHICH  IS  RELATED  TO  THE  DEFLECTION  AT  THE 

C      END  .  THE  SLOPE  AT  THE  MIDDLE  ,  WHICH 

C      IS  ZERO  BY  5YMHETRY,MUST  BE  IMPOSED. 

C      FOR  OTHER  TYPES  Cif    BOUNDARY 

C      CONDITIONS, THIS  SUBROUTINE  IS  NOT  TRUE. 


C 


IF{-N.EQ,U  GO  TO  3 

KK=IDBC(N) 

DO  4  K=l,NNOD 

SYSK(KK,K>  =0,CD<-00 

SYSK{K,KK)=0. 00+00 

DO  4  I  =  l,ivJNuD 

DO  4  J=1,NN0D 

SYSC{KK, I , J) ^O.OD+00 

SYSC(KK,J,  n=O.OD*-00 

SYSC(  KK,!<K,KK)=O.OD+00 

SYSF(K.K)=O.OD  +  00 

SYSK(KK,KKi=l.D+00 

CONTINUE 

RETURN 

END 
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SUBROUTINE  CURF I T{ LP , PO,B) 
C 

C      USING  LSQPL?  LIBRARY  SUBROUTINE  ,  THE 
C      IRRATIONAL  J^OWER  CURVE  WITH  P  NOT  EQUAL  TO 
C      AN  INTEGER  IS  REPLACED  BY  A  LEAST  SQUARE 
C      BEST  FIT  SECOND  ORDER  POLYNOMIAL 


C 


IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  8( 10) ,WI (100) ,Y( 100) ,DELY( 100) f 
1      POdO;  tSB(100),F2(  iOO)  ,X(  100)  ,ER(100) 
REALMS  TITLE(10)/10-'     •/ 
DATA  Ml/iOO'n.DO/ 
NMAX=-2 
P=PO(LP) 

IF(P.EQ.2.D+00)  GO  TO  10 
N  =  5 

DO  2  1=1, N 
FLOATI=I 
X(I>=FLGATI 
2  F2(I)=X(  D^^P 

CALL  LSGPL2(N,iNMAX,X,F2tWI,Y,DELY,B,SB,TITLE) 

DO  6  1=1, N 

IF(F2(  n.LE.l.D-06)  GO-TO  7 

ER(I)  =  DELY(I)/F2(I  ) 

GO  TO  6 

7  ER(I)=O.OD+00 
6  CONTINUE 

£RR0R=ER(1) 

DO  8  1=2, N 

IF(DAaS(ER(I)) .L£. ERROR)  GO  TO  8 

ERROR  =  DABS(ER(  I)  ) 

8  CONTINUE 
PRINT  9, ERROR 

9  FORMAT(//, lOX, • THE  MAX  RELATIVE  ERROR  IS=*,015.7) 
GO  TiD  11 

10  B(l)=0.0D<-00 
B(2)=0. 00*00 
B(3)=l.D+00 

11  RETURN 
END 
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SLBROUTINE  FSOIL(B) 

IMPLICIT  REAL*8( A-H,0-Z) 
C 

C      THE  EXPERIMENTAL  DATA  ON  A  TYPICAL 
C      NATURAL  SO  I L ,  I . E  .,  RE  ACT  I  ON  OF  FOUNDATION 

C      VERSUS  LOADING  CONDITION  FROM  TABLE 

C      AND  FIGURE ,ARE  APPRIMATED  BY  A  LEAST 

C      SQUARE  BEST  FIT  SECOND  ORDER  POLYNOMIAL 

C      USING  LSQPL2  LIBRARY  SUBROUTINE 

C 


c 


1 


DIMENSION  B( 10) ,WI( 10) tY(lO) ,DELY(10) , 

SB( 10),F2( 10) ,X( 10J,ER(10) 
REAL*8  TITLE( 10)/10-»     •/ 
DATA  WI/IO^I.DO/ 


REAL  V.  ..,w^..^. 
DATA  WI/lO^l.DO/ 
NMAX=-2 
N  =  6 

DO  6    1=1, N 
2  READ(5,3)  X( I) tF2( I) 


3 


READ(b,3)  XlIJtFZd) 

FORMAT( 2G13.6) 

CALL  LSUPL2(N,NMAX,XiF2,WI,Y,DELY,B,SB,TITLE) 


C 

C  THE  SEA  SEDIMENT  FOUNDATION  IS  MODELED 

C  FROM  INLAND  SOIL  DATA  BY  A  SCALE  1:3000 

C  OMIT  THE  Three  following  cards  IN  CASE 

C  IN  CASE  OF  INLAND  FOUNDATION 


B(U=6i  1)*  8.  30+0  0/30  00. 0+00 
B  (  2)  =B  ( «^ )  =^8,  3D  f  00/3000.0  +  00 
6(3)  =  B(  3  )=<=o.  3D +00/ 3000.0  +  00 
RETURN 
END 
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SUBROUTINE  VGUES  S  (  NNOD,  NEL  EM,  V  SCALE,  IT,MAX,V) 
C 

C      THIS  SUBROUTINE  PROVIDES  AN  INITIAL 
C      ESTIMATE  FOR  NONLINEAR  PROBLEMS 
C 

IMPLICIT  REAL*3( A-H,0-Z) 

DIMENSION  V(20) 

ITMAX=100 
C 

C      THE  MAIN  PROGRAM  IS  DESIGNED  IN 
C      SUCH  A  WAY  THAT  IF  ONE  INITIAL  ESTIMATE 
C      DOES  NOT  GIVE  SOLUTIONtTHE  SCALE  WILL  BE 
C      CHANGED, AND  THE  ITERATION  RESTART  AGAIN. 
C      HOWEVER, AFTER  10  CYCLES, IF  NO  SOLUTION  IS 
C      OBTAINED, THE  READER  MAY  CHANGE  THF.  INITIAL 
C      SCALE  BY  HAND  IN  THE  SUBROUTINE  INCHK 
C 

VSCAL E=l. 50+00* V SCALE 

TAU=O.0u+OO 

NT=NELEM+l 

DC  22  1=1, NT 

11=2*1 

111=2*1-1 
C 

C      FOR  SIMPLICITY  ,  A  SECOND  ORDER  POLINOMIAL 
C      IS  USED  FOR  ESTIMATE  SOLUTION  .  THE  READER 
C      MAY  DESIGN  OTHER  ESTIMATE  FUNCTION 
C      AS  HE  WHISHES 


C 


V(II )=VSCALE*2.D+00*(4.D+00-TAU) 
V(III)=VSCrtLtv(  (  8.D  +  00*TAU-TAU**2) ♦a.OlD+00) 

22  TAU=TaU+1 .0+00 
Vil)=0.0D+00 

PRINT  2i, (V( I) , 1=1 ,NNOD} 

23  FORMAT( ////,15X, • INITIAL  GUESS  : « ,/// » d 5X , 01 5 .6 , / ) ) 
RETURN 

END 
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SUBROUTINE    DIFSOL 
C 

C  USING    FINITE    DIFFERENCE    METHCO    TO    SOLVE 

C  NCMINEaR    PROBLEMS, THIS    SUBROUTINE    IS 

C  USED    STRICTLY    FOR    THE    CASE    OF    UNIFORM    BEAMS 

C  WITH    UNIFORM    LOADS    AND    THE     REACTION    OF 

C  FOUNDATION    IS    PROPORTIONAL    TO    THE       SQUARE 

C  OF    ITS    DEFLECTION 

C 

IMPLICIT    REAL*8( A-H,0~Z) 

COMMON    EK(4t4),DC(4,4),DM(4,4),DN(4,4,^), 
i  ICORR(20,4)  .SYSK( 23, ^0)  ,SYSC( 20,20,20)  , 

2  SYSF(20)  ,EFU),IDBC(20)  ,8(10), 

3  NNOO,NELEM,NDOF,NBC,TL,AAl, AA2, 

4  X, EI ,Q, Alpha, SS0IL,A1,A^ 
DIMENSION    U(I0),WD(100) 

c 

C      TWO  BEAM  MODELS  ARE  CONSIDERED  : 

C      EXTERNAL  DIFF4  CORRESPONDS  TO  4-ELEMENT 

C      MODEL, AND  EXTERNAL  0IFF8  TO  8-ELEMENT  MODEL 

C 

EXTERNAL  DIFF4 

EXTERNAL  DIFF8 
1  FORMAT ( • 1*  ) 

ND=NELEM 

OD  =  ND 

NSIG=5 

EPS=l.D-06 

NCCUNT=1 

USCALE=1000.D<-00 

3  ITMAX=100 
NC0UNT=NC0UNT+1 
TAU=O.0D+0O 
USCALE=1.5D<-00*USCALE 
DO   4    1=1, NO 
TAU=TAU+2.DfOO/DD 

4  U( I)=USCALE*(TAU'.25D+00*TAU**2) 
IF(ND.GT.2)    GO    TO    5 

CALL    ZSYSTM(DIFF4,EPS,NSIG,ND,U, ITMAX , WD ,PAR, I ER) 
GO    TO    8 

5  CALL    i:SYSTM(DIFF8  , EPS, NSIG,ND,U,  ITMAX, WD, PAR,  lER) 

8  IF(NCOUNT.GT.30)     GO    TO    6 
IF(IER^NE.O)     GU    TO    3 

6  PRINT    7, lER, ITMAX 

7  F0RMAT(////,i5X, » I ER= • , I 5//15X , • ITMAX=« ,15) 
DO    9    i=I,NU 

9  U( I )=U( i)/EI 
PRINT  10 

10  FORMATS ////, 15X, 'F.D.M.  RESULT*,////) 
PRINT  1 1, (U(  I)  ,1=1  ,ND) 

11  F0RMAT{15X,G20.6/) 
RETURN 

END 
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SUBROUTINE  STIFFl 
C 

C      THIS  SUBROUTINE  IS  USED  TO  COMPUTE  THE  ELEMENT 
C      STIFFNESS  MATRICES  AND  THE  SYSTEM  STIFFNESS 
C      MATRIX  FOR  LINEAR  PROBLEMS 
C 

CCMMON     EK(^t4) ,DC(4,4) ,DM( >,4} f DN( A,4,4I » 

1  IC0RR(20,4)  ,SYSK( 20,20)  ,SYSC(20»20,2C) , 

2  SYSF(20),EF(4)  ,  IDciC(aO)  ,B(IO)  , 

3  NNOD,  NELEMffJUOF  tNBC  ,  TL  ,  AAl ,  AA2, 

4  X,EI,QTALPHAtSSOlL,AI, A2 
ALPHA=ALPHA/£I 
SSCIL=SS0IL/EI 

EK(i, I)=12.*(I./X**3) 

EK(l,2J=6.-( l./X-*2) 

EK(l»3)=-i2.*(l./X*=i=3) 

EK(1,4)=6.*(  l./X=«==^2) 

EK(2,2)=4.-(1./X) 

EK(2,3)=-6.^=(1./X**2) 

EK(2,4J=2.-( l./X) 

EK(3,5)=12.-{1./X**3) 

EK(3,4)=-6.-( i./X**2) 

EK(4,4)=4. =*=(!. /X) 

EK(2,I )=EK(l,2i 

EK{3,l)=EK{i,3) 

EK(3,2)=  EK(2,3) 

EK(4,1) =EK(1,4) 

EK(4,2}=EK(2,4) 

EK(4i3i=EK(3,4) 
C 

C      THE  FOLLOWING  ARE  THE  ELEMENT  STIFFNESS 
C      MATRICES  ASSOCIATED  TO  THE  FOUNDATION 


C 


DC(l,li=1.2/X 
DC(1,2J=0.1 
DC(1,3)=-L)C(  If  1) 
DC{l,4)=DC(I,2i 
DC(2,1)=DC(I»2) 
DC(2,2)  =0.133333=^X 
DC(2,3)=-DC( 1,2) 
DC(2,4)=-0.G333:»3*X 
DC(3,1)=DC(1,3) 
DC(3,2)=DC(2,3) 
DC(3,3)=DC(I,1 ) 
DC(3,4)=DC(2,3) 
DC(4,1)=0C(1,4) 
DC(4,2)  =[)C(2,4) 
DC(4,3)=DC(3,4) 
DC(4,4) =DC(2,2) 
0M(  It  1)  =  .3  71429=:--X 
DM(l,2)  =  .a52381=^X*«2 
DM(i,3)=.12d571-X 
DM(1,4)=  -.030952'!=X**2 
DM(2,2)=.009524*X**3 
DM(2,3i=  .O30952*X*'?2 
0M(2,4)=-.007I43>i'X**3 
DM(3,3)=.37i^29-X 
DM(3,4)=  -.05238i*X**2 
0M(4,4)=  .009524*X=<'*3 
DM  I  2, 1 )=DM( 1,2) 
DM(3,1>=  DM( 1,3) 
DM{4,i) =  DM( 1,4) 
DM{3,2)=  OM(2,3) 
DM(4,2)=D,^(2,<t) 
DM(4,3)=  0M(3,4i 
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c 
c 
c 
c 
c 


c 
c 
c 
c 
c 


SINCE  ALL  ELEMENT  STIFFNESS  MATRICES  ARE 
ASSOCIATED  TO  THh  LINEAR  TERMS  ,  THEY 
MAY  BE  ADOED  ^S  FOLLOWS 

DO  1   1=  1,ND0F 
DO  1   J=  1  ,NDOF 

1  EK(1  ,  JJ=EK(  i  ,  J)+ALPHA*DM(  I  ,  J  )  f  SSGI  L'i'DC  (  I 
DO  2  I=1,NNGD 

DO  2  J  =I,NNOD 

2  SYSK( I , J>=  0.0 


4 
3 


J) 


THE  SYSTEM  STIFFNESS  MATRIX  IS  OBTAINED   BY 
COLLECTING  THE  STIFFNESS  MATRICES  OVER  ALL 
ELEMENTS   OF   THE   SYSTEM 


DO  3 
DC  4 
NS1  = 
DO  4 
NS2= 


N=     ItNELEM 
i=l,NDQF 
ICORR{N, I) 
J=     ItNOOF 
ICCRR(N, J) 

SYSK( NSi,NS2)= 

CONTINUE 

RETURN 

END 


SYSK(NS1,NS2)    ■•■    EK(I,J) 
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SLBROUTINE  SOLVE ( SYSK, SYSF , V»NNOD ) 
C 

C-      BASED  UN  GAUSSIAN  ELIMINATION  TECHNIQUE 
C      THIS  EQUATION  SOLVER  IS  USED  TO  SOLVE 
C      THE  LINEAR  ALGEBRAIC  SYSTEM  ONLY 
C 

DIMENSION  SYSF(20)  ,SYSK(20,20) , 
1  IV(2D» ,SYSFPi20) ,V(ZO) 

DC  I  K=1,NN0D 
1  IV(K)=K 
C 

C      GAUSSIAN  ELIMINATION 
C 

NMi=NNOD-l 
DC  1004  1=1, NMl 
IP1=UI 

ATCP=AbS(SYSK( 1,1)) 
IPkIME=I 

DO    1100     IP=IP1,NN0D 

IF(A5S( SYSK(  IP, I  ) i  .LE.ATGPi    GO    TO    1100 
AT0P=A8S<SYSK( IP,1 )) 
IPRIME=IP 
1100    CONTINUE 

IF( IPRIME.EO.I )  GO  TO  1200 
ISAV=IV( IPklME) 
IV(IPRIME}=IV( I) 
IV(I }=ISAV 
DO  1300  J=ltNNGO 
TEMP=SYSK( IPRIME, J) 
SYSK( IPRIME, J)=SYSK(I , J) 
1300  SYSK( I, J)=TEMP 
1200  DO  1005  K=IP1,NN0D 

FACTG}<=SY5K(K,I  )  /SYSK(  I,I> 
SYSK(K, I)=FACTDR 
DO    1006    J=IP1,NN0D 

SYSK(K, J)=SYSK(K, J)-FACTOR*SYSK{ I, J) 
1006    CONTINUE 
1005    CONTINUE 
1004    CONTINUE 

DO    9    I=1,NN00 
K=IV(n 
9    SYSFPd  )  =  SYSF(K) 
DC    10    J=1,NM1 
JPl=Jf 1 

DO    11    K=JP1,NN0D 

SYSFP(Ki=SY$FP(K)-SYSK(KT J)*SYSFP(J) 
11    CONTINUE 
10    CONTINUE 
C 
C  BACK  SOLUTION 


C 


SYSF(NNCD)^SYSFP(NNGD)/SYSK(NNOD,NNODJ- 

DO    13    II=i»NMl 

I=NNOD-Ii 

IPl=H'l 

SUN^O.O 

DC  14  J=IP1tNN0D 
14  SLM=SUM^-SYSK(  I  t  J  )*SYSF(  J) 
13  V(l)=(SYSFP( I )-SUMi/SYSK( I , I) 

RETURN 

END 
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SUBROUTINE  RESULT ( NNOD ,SY SF, £1 , I ER, ITMAX , NCODE ) 
C 

C      THIS  SUekCUTINE  IS  USED  TO  PRINT  THE 
C      SGLUTICN  UF  EITHER  LINEAR  OR  NGNLINEAR 
C      PRObLEMS.FGR  LINEAR  PROBLEMS  ,  OMIT 
C      'IMPLICIT  HEAL<^8(  A-Ht  C-Z)  •  CARD 
C 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  SYSF(20) ,V(20) 

IFCNCOOc.EC.l J  GO  TO  19 

PRINT  28  ,  lER, ITMAX 
28  F0B,MAT(//,15X,«  IER  =  '  ,1X,I2,////, 
ili)X,MTMAX  =•  ,iX,I3»////) 

19  NN=NNuD/2 

DC  20  I=lfNN 
11  =  2*^1 
JJ=2=<^i-'l 

20  V(I)=-SYSF(  JJI 
PRINT  3  0 

30  FCR^iAT(////,  15X,  •DEFLECTION  :',//) 
PRINT  31,( V( li ,i=l .NN) 

31  FORMAT (  5X,G20.6,/) 
C 

C      IF  THE  DEFLECTIONS  AS  WELL  AS  THE  SLOPES 

C      AT  THE  NUDAL  POINTS  ARE  DESIRED  f 

C      ADD  THE  FOLLOWING  CARD 

C      PRINT  31,(SYSF(  U  ,  I=^l,NNOD) 


C 


RETURN 
END 
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C  THE    FUNCTION    AUX    WHICH    IS    TRANSLATED 

C  FROM    EQUATIOiM PAGE ,GLi^EftALLY 

C  REMAINS  THE    SAME. BUT    THE    READER    HAS    TO 

C  REDESIGN    OTHER    EXTERNALS    t     SUCH    AS    QFUNCT 

C  EFUNCT,SFUNCT,bZ,63fWHICH    MAY    VARY    , 

C  DEPENDING    CN    LOADING    COND I T IGNS , BE AM 

C  CHARACTERISTICS    AND    TYPES    OF    FOUNDATIONS 


C 


DOUBLE    PRECISION    FUNCTION    AUX{V,K,PAR) 

IMPLICIT    REAL*8( A-H,0-Z) 

CCMMCi\    £K(4,^)  ,DC(4,4)»DM(4t4)  »DN(4t4,4)  , 

1  IC0kR(20,4I  » SY ski  20,20) ,SYSC( 20,20, 20) t 

2  SYSFt  2G),EF (4), I0BC(20) ,6( 10), 

3  NNOD,NELEM,NUOF ,NBC , TL , AA 1 , AA2, 

4  X, EI, g, ALPHA, SS0IL,A1  ,A2 
DIMENSION    V(20) 

SKT£MP  =  0.0D<-00 
SCTEMP=0.0D+00 
DO    1    1=1,NNCD 

SKTEMP=$KTEf^.P+SYSK(K,I  i  ^\J  {  I) 
DO    I     J=l,NNOO 
1    SCTEMP=SCTEMP+SYSC(K, I , J ) * V ( I) *V(  J) 
AUX=SKTE.MP+SCTEMP-SYSF{K) 
RETURiM 
END 

DOUBLE  PRECISION  FUNCTION  QFUNCT(N,X) 

IMPLICIT    REAL*8( A-H,0-Z) 

FLOATN=N 

XI  =  (FLCATN-.5D  +  D0)=!=X 

IF(N,NE.l)    GO    TO    3 

QFUNCT  =  0.0D<-00 

GO    TO    4 

3  QFUNCT=(  {3.  50  +  00/3. D+00)  +  (  I  .0  +  0  0/ 150.D4-00)  * 
1  XI)=:-l.p  +  04 

4  RETURN 
END 

DOUBLE  PRECISION  FUNCTION  B2(NtX) 

IMPLICIT  REAL*8( A-H,0-Z) 

FLOATN^N 

XI  =  (FLGATN'.5D  +  00i'i-'X 

IF(N.NE.l)  GO  TO  3 

B2=2000.0+0G+40.D*00*XI 

GO  TO  4 

3  62=(  10G-0i-0U/l5.D4-00J«XI  +  (  1 1000.  D+00/3  .D  +  OOJ 

4  RETURN 
END 

DCUdLE  PRECISION  FUNCTION  B3(N,X) 

IMPLICIT  REAL-31 A-H,0-Z) 

FLOATN=N 

Xl=( FLOAT N-.5D+00)*X 

B3  =  .  50  +  00^^X1+300.0 +00 

RETURN 

END 

DOUBLE  PRECISION  FUNCTION  EFUNCT(N,X} 

IMPLICIT    REAL=^8(  A-H,0-Z) 

FL0ATN=N 

X1=(PL0ATN-.5D+00)*X 

IF{N.Nc. I)     GO    TO    3 

EFUNCT=i-D+li+( I .D+11/50.D+00) «XI 

60    TO    4 

3  EFUNCT=2.0+ll 

4  RETURN 
END 

DOUBLF  PRECISION  FUNCTION  SFUNCT(N,X) 

IMPLICIT  REaL*8( A-H,0-Z} 

SFUNCT=0. 00+00 

RETURN 

END 
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DOUBLE  PRECISION  FUNCTION  DIFF4( Ut K, PAR) 

IMPLICIT  RL-AL=5=8(  A-H,0-Z) 

CCMMGN  EK(4,^)tDC(4,4),0M(^,4),DN(4,4,4)t 

1  I  CORK <^ 0,4)  ,SYSK( 20,20) ,SYSC( 20,20,20) , 

2  SySF( 20) ,EF (4) ,IDBCl20i ,B( 10), 

3  NNCD,NELEM, ND OF ,N8C , XL , AAI , AA2, 

4  X, EI, Q, ALPHA, SSGIL,Ai ,A2 
DIMENSION  UCIO) 

GO  TO  (5,10) ,K 
5  DIFf  4=  1,53 6L' +00 *U(  1 ) +0.  OOID  +  OO*  A2*(TL**4  )*U  ( 1  )**2- 
11.02^D+00*U(2)-O.OOID+00-Q-(TL**4J 
RETURN 
10  DIFF4=-2.048D<-00*U(1)  +  1.536D+00*U(  2)  + 

10.0GlD<-00*A2-{  TL--4)*U{2)*-2-0.  00iD*-00*Q*(  TL**4) 
RETURN 
END 


DOUBLE 

IMPLICI 

CCMMCN 

1  IC 

2  SY 

3  NN 

4  X, 
DIME  MSI 
GO  TO  ( 

5  D1FF8=5 

1  (X*-4)  <' 

2  -Q* 
RETURN 

0  DIFF8=- 
1U(4)+(X 
2     -Q- 
RETURN 
15  DIFFS= 
1(X**4)=^ 
2     -Q* 
RETURN 
20  DIFF8= 
1 (X«*4)* 
2     -Q« 
RETURN 
END 


I 


PRECISION  FUNCTION  D IFF8 ( U,K, PAR ) 

T  REAL*8( A-H,0-Z) 

EK(4,4),0C(4,4),DM(4,4),0N(4,4,4), 

URK(20,4) , SYSK( 20,20 ) ,SYSC( 20,20, 20) , 

SF(20) ,EF{4)  ,IDBC(20) ,B(10) , 

OD, NELEM, NDOF ,NBC , TL , AAI , AA2  , 

EI, Q, ALPHA, SS0IL,A1,A2 

ON    U( 10) 

5,10, 13,2 0),K 

.0+GO^Ui  1  )-4.D+00*U(  2)     +    U(3)     4- 

(  A1*B(  l)4-AAi-U(  1)<-AA2=«'U{  1J**2) 

X**4 

4.D*00=s=U{  1  )+6.D<-00*U(2)-4.D  +  00=^U(  3  )  <■ 
**4)*(  Al^-3(  l)<-AAl*U(2J+AA2-U(2»-*2) 
X**4 

U(  1  )-4.D  +  00=i^U(2  )+7.0  +  00*Ul3)-4.D<-00*U(4)4- 

(  Al<^B(  l)  +  AAl«U(3)+AA2*U(3)*-2) 

X=i'*4 

2.D+00*U(  2)-8.Df00*U(3)+6.D  +  00*U(4)<- 
(  Al^f^B  (  1 )  <■  AA  1*U  (  4  ) +AA2^U  (  4  )**=2  ) 
X=5'=<=4 
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